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ABSTRACT 
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The  development  of  a  numerical  procedure  for  the  treat- 
ment of  nonsimilar,  unsteady,  laminar  boundary  layers  is 
presented.   The  solution  is  obtained  from  the  laminar,  iso- 
thermal, incompressible  boundary  layer  equations  employing 
a  modification  of  the  integral  matrix  procedure  of  Bartlett 
and  Kendall.   Solutions  of  example  problems  are  presented 
for  steady  Blasius  and  Howarth  flows,  and  for  oscillating 
Blasius  flow.   Agreement  with  the  known  classical  results  is 
satisfactory  and  establishes  the  general  feasibility  of  the 
method.   Core  storage  requirements  of  130,000  bytes  allow 
consideration  of  as  many  as  25  nodal  points  across  the 
boundary  layer,  50  time   increments  per  oscillation  cycle 
and  50  streamwise  stations.   Solution  of  oscillating  Blasius 
flow  considering  8  nodal  points  and  16  time  increments 
requires  13.49  seconds  for  one  streamwise  station  utilizing 
IBM  360/67  time  sharing  capabilities. 
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I.   INTRODUCTION 

A.  GENERAL 

All  aerodynamic  flows  are  to  some  degree  influenced  by 
unsteadiness.   Rotating  stall  in  turbomachinery ,  aerodynamic 
stall  flutter  and  hydrofoil  flow  fields  are  but  a  few 
practical  problems  in  which  the  nonsteady  fluid  flow  phe- 
nomena are  of  prime  importance.   Additionally,  the  recent 
emphasis  on  the  development  of  rotary  wing  aircraft  has 
created  a  new  demand  for  knowledge  of  the  behavior  of 
boundary  layers  in  unsteady  flow. 

B.  ANALYTICAL  HISTORY 

Early  analytical  treatments  by  Stokes,  Rayleigh  and 
Schlichting  (Ref.  1)  considered  only  a  special  case  of  the 
general  problem:  that  of  unsteady  viscous  fluid  flow  with 
no  mean  flow  or  pressure  gradient.   Lighthill  (Ref.  2)  made 
the  first  significant  investigation  of  small  sinusoidal 
oscillations  superimposed  on  a  mean  flow.   Lighthill' s 
analysis  was  based  on  a  small  perturbation  treatment  of  the 
boundary  layer  equations  for  sinusoidal  flow,  with  only 
the  first  order  terms  retained.   The  determination  of  the 
existance  of  the  "quasi-steady"  regime  at  low  frequencies, 
and  the  "shear  wave"  regime  at  higher  frequencies  was  the 
principal  result  of  this  analysis.   Glauert  (Ref.  3)  extended 
Lighthill' s  work  by  considering  the  boundary  layer  in  the 
vicinity  of  a  stagnation  point.   Carrier  and  DiPrima  (Ref.  4), 
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using  a  system  of  equations  derived  using  a  modified  Oseen 
linearization,  confirmed  the  existance  of  the  "shear  wave" 
solution  and  the  phase  advance  of  the  wall  shear  determined 
by  Lighthill.   Nickerson  (Ref.  5)  initially  retained  per- 
turbation terms  beyond  the  first  order,  thus  widening  the 
scope  of  Lighthill' s  work. 

In  addition  to  his  analytical  development,  Nickerson 
made  hot-wire  measurements  of  the  laminar  boundary  layer  on 
a  flat  plate  oscillating  in  a  Blasius  mean  flow.   His  work 
partially  confirmed  both  his  own  and  Lighthill' s  analyses, 
but  mechanical  difficulties  precluded  the  completion  of  the 
experimental  program.   Hill  and  Stenning  (Ref.  6)  imposed 
sinusoidal  oscillations  directly  on  mean  Blasius  and  Howarth 
flows  and  were  able  to  confirm  the  analytical  results  of 
both  Lighthill  and  Nickerson.   Recently  an  integral  analysis 
of  the  unsteady  boundary  layer  was  introduced  by  Koob  and 
Abbott  (Ref.  7) ,  who  succeeded  in  demonstrating  the  method 
for  the  case  of  an  accelerating  flat  plate. 

C.   APPLICATION  OF  NUMERICAL  METHODS 

With  the  advent  of  the  high-speed  digital  computer,  it 
has  become  possible  to  apply  numerical  methods  to  the 
analytical  studies  of  viscous  fluid  flows.   Fromm  (Ref.  8) 
and  Chorin  (Ref.  9)  introduced  modified  finite  difference 
methods  for  the  treatment  of  viscous  flows  which  could 
feasibly  be  applied  to  the  boundary  layer  problem.   Smith 
and  Clutter  (Ref.  10)  and  Kleinstein  and  Nabi  (Refs.  11  and 
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12)  applied  finite  difference  methods  to  the  steady  boundary 
layer  equations.   Successive  approximation,  linearization 
and  integral  methods  have  also  been  applied  to  the  boundary 
layer  equations  (Refs.  13,  14  and  15).   A  major  drawback  of 
most  of  the  methods  mentioned  is  that  in  order  to  reach  a 
solution  to  a  given  problem  large  computing  times  are  nec- 
essary. 

The  introduction  of  the  Boundary  Layer  Integral  Matrix 
Procedure  by  Bartlett  and  Kendall  (Ref.  16)  provided  a  major 
breakthrough  in  the  field  of  fast  computer  solution  of  the 
steady  flow  boundary  layer  problem.   The  method  was  developed 
as  a  means  of  obtaining  rapid  computations  in  the  laminar 
regime  for  steady,  multicomponent ,  chemically  reacting 
boundary  layers;  and  the  logical  extension  of  the  technique 
was  to  encompass  the  unsteady  flow  condition. 

The  purpose  of  this  investigation  is  to  develop  a  com- 
puter program  capable  of  solving  the  unsteady,  laminar 
boundary  layer  problem  by  adapting  the  rapid  solution  tech- 
niques of  the  Boundary  Layer  Integral  Matrix  Procedure. 
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II,   ANALYSIS 

A.   GOVERNING  EQUATIONS 

In  this  development,  utilizing  the  integral  matrix 
method,  it  is  convenient  to  make  the  assumption  of  incom- 
pressible isothermal  flow.   The  resulting  laminar  boundary 
layer  equations  for  two-dimensional  (n=0)  or  axisymmetric 
(n=l)  flow  are  as  follows: 

Continuity: 

i   8 (r  nu)    . 

-L-   °   +  |x=  o  (i) 

r  n     ds  9y  v  ' 

o 


Momentum 

9t  + 

9u 
U9? 

+ 

9u 

va7 

1 

P 

9P 
9  s 

+ 

U 

P 

9  u 
9  2 

y 

H  + 
9t 

U9t 

= 

1 

P 

9P 

9s 

(2) 


(3) 


where  s  and  y  are  respectively  the  streamwise  and  normal 
components,  u  and  v  are  the  corresponding  velocity  compo- 
nents, U  is  the  streamwise  velocity  component  outside  the 
boundary  layer,  r   is  the  radius  of  the  body  perpendicular 
to  the  axis  of  revolution  for  an  axisymmetric  body,  n  is 
zero  for  a  two-dimensional  body  and  unity  for  an  axisymmetric 
body,  t  is  time,  p  is  density,  y  is  viscosity  and  P  is 
pressure . 
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B.   BOUNDARY  CONDITIONS 

The  boundary  conditions  considered  for  the  governing 
equations  as  stated  are: 


u      = 

w 

0 

V        = 

w 

0 

u       = 
e 

U(s 

,t) 

8u 
e 

9y 

=    0 

(4) 
(5) 
(6) 

(7) 


where  the  subscripts  w  and  e  refer  respectively  to  wall  and 

boundary  layer  edge  conditions.   Equation  (4)  is  the  require- 
ment of  zero  slip  at  the  wall,  equation  (5)  implies  no  mass 
transfer  at  the  wall,  and  equations  (6)  and  (7)  are  repre- 
sentative of  the  assumption  that  the  boundary  layer  edge 
velocity  is  equal  to  the  free  stream  velocity. 
The  free  stream  velocity  may  be  written: 


U  =  (U   -  b,s)  (1  +  e  cos  u)t)  (8) 


where  U   is  the  mean  free  stream  velocity  at  the  leading 
edge,  b,  is  the  rate  of  change  of  the  mean  free  stream 
velocity  with  distance  from  the  leading  edge,  z    is  the  ratio 
of  the  amplitude  of  the  free  stream  oscillation  to  the  local 
mean  free  stream  velocity  and  w  is  the  frequency  of  the 
oscillation. 
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C.   TRANSFORMATION  OF  EQUATIONS 

Paralleling  the  development  of  Bartlett  and  Kendall 
(Ref.  16)  a  transformation  of  variables  is  performed  util- 
izing a  combination  of  the  Levy  and  Mangier  and  the  Howarth- 
Dorodnitsyn  transformations  which  is  known  by  a  variety  of 
names  including  Levy-Lees,  Lees-Dorodnitsyn  and  Mangler- 
Dorodnitsyn.   This  transformation  is  as  follows: 


5  -  /   uepeyero2nds  (9) 

o 


r  nu   y 
n  =  — — -  /  P  dy  (10) 


/2?   o 

Simplifying  equations  (9)  and  (10)  to  encompass  the 
assumptions  of  incompressibility  and  constant  temperature 
the  transformations  become: 


5  =  py  /   uerQ2n  ds  (11) 

o 


n 
r   pu 

-2 ^y  (12) 


•25 

The  transformation  to  Levy-Lees  variables  is  greatly 
implified  by  use  of  the  operator: 


■Ni1  ♦  ^  -  «Vo°  ■"*  «ty -$+$*#■>    (13> 


where  f  is  the  stream  function  defined  as: 


f  -  f«  ■  "  £■  dn  (14) 

o   e 


the  prime  represents  partial  differentiation  with  respect  to 
n  so  that: 


(15) 


and   a*   is    a  normalizing  parameter: 


Pliur" 
a*   =   (16) 


/2C 

In  addition,  the  partial  derivatives  are: 

3(  )    a* 


3y    y 


(  )  (18) 


^— v-  =  (V)    (  }  (19) 


Applying  equation  (15)  to  the  velocity  in  the  partial 
derivative  with  respect  to  time  yields: 

3u       3f     -;i    e  /OA> 

at  "  ue  "St  +  f  TF  (20) 

Using  equations  (13)  through  (20)  and  applying  equation 
(6)  to  equation  (3)  results  in  the  following  transformed 
momentum  equation : 


Ill       II  I  •£ 

f  +  ff  +  3(1  -  f   ) 


•2  If'  -iil  -f"  af  )-  -^-  f9f'  +  ff'-i)3AnUe>  -o  ran 


where 


f 

w 

" 

0 

f 

w 

= 

0 

1 

f 

e 

= 

1 

dinu 

-  2  -sih|  (22) 


Transforming  the  boundary  conditions,  equations  (4) 
through  (7)  become: 

(23) 
(24) 
(25) 
fe  =  0  (26) 

D.   TAYLOR  SERIES  EXPANSIONS 

At  any  streamwise  station  s,  the  boundary  layer  may  be 
divided  into  N-l  strips  which  are  joined  by  N  nodal  points 
n • ,  where  i  is  1  at  the  wall  and  N  at  the  boundary  layer 
edge.   The  primary  dependent  variables  f.  and  their  normal 

derivatives  are  related  by  Taylor  series  expansions  such 

i 

that  the  f.  are  represented  by  connected  cubics  with  con- 
tinuous first  and  second  derivatives  at  the  nodal  points. 
This  is  commonly  called  a  spline  fit. 

Expanding  the  primary  variables  in  Taylor  series  and 

i  H 

truncating  at  f.   the  following  linear  equations  result: 
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.f.+l+fi+f:a,  +  f:^+<"in!+£:;i^r=°    <27» 


i  i  u  ,MfSn^  MI        fin2 

fi+1  +   f ±  +   f .    «n  +   f i  %-  +  fi+1  22_  =  o  (28) 


ii  ii  III       r  ill  £„ 

fi+i  +  £i  +  fi    T  +  fi+i  f  "  °  (29) 


«n  =  ni+1  -  r»t  (30) 


III  Ml 

-  m        fi+l        -        ^  (31) 


Equations     (27)    through    (29)    represent    a   set   of    3(N-1) 

i 

equations  where  f  is  represented  as  a  quartic,  f   as  a  cubic, 

f   as  a  quadratic  and  f    is  considered  to  vary  linearly 
between  n.  and  n  •_.-[• 

E.   AXIAL  AND  TIME  DERIVATIVES 

The  axial  derivatives  of  equation  (21)  are  treated  as 
logarithmic  two  or  three-point  differences  such  that: 

2  \mi\   -  do(  h  +  V  >i-i  +  <V  h-2  (32) 

where  for  two-point  differences: 


d^  =  — r2- d,  =  -d^         d   =  0         (33) 

o    £A^  1     o  2 
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and  for  three-point  differences: 


,  _  „£A£-l+£A£-2         9 
d  =  2. — r r a,  =  -2 — r 

o    „  A„  .,  „  A„  ^     1     n  A 


£-2 


£-1  £  Jt  — 2 


£  £-1  £-1  £-2 


d2=2^ 


£A£-1 


£  £-2  £-l"£-2 
(34) 
In  equations  (32)  through  (34),  (  )  .  ,implies  the  previous 
streamwise  station. 


£A£-1  "  inh    ~    "nh-l    =  *n(V^£-l> 


(35) 


Differencing  the  time  derivatives  in  the  same  way: 


flLi 

L  dt 


■  To(  >k  +  Tl(  >k-l  +  T2(  >k-2 


(36) 


where  for  two-point  differences: 


o    n  A 


k"k-l 


T  _  _T 
1     o 


T2  =  0 


(37) 


and  for  three-point  differences: 


T  _  kAk-l+kAk-2      _   _ 
o   ,  A,  ,  ,  A,  n         ll        ,  A 


kAk-2 


kAk-l 


-1   k-lAk-2  2      kAk-2    k-lAk-2 


(38) 


In  equations  (36)  through  (38),  (  ),  ^implies  the  previous 
time . 


kAk-l  =  fck  "  Vl 


(39) 
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Equations  (37)  through  (39)  may  be  greatly  simplified 
by  considering  the  utilization  of  equal  time  spacings.   The 
resulting  equations  for  the  treatment  of  time  derivatives 
over  equal  time  spacings  are  as  follows: 

Two  point  differences: 

To  "  6T  Tl  =  "To  T2  ■  °         (40) 

Three-point  differences: 

T   =   1  T   =  — —  T   =  — —        (41) 

o    26t  1     St  l2        26t        K-      } 

where : 


St   m   fck  "  **-l  "   Vl  "  \-2  (42» 


F.   INTEGRATION  OF  THE  BOUNDARY  LAYER  EQUATION 

The  transformed  streamwise  momentum  equation,  equation 
(21),  is  integrated  at  constant  £  from  n-_-,  to  n .  ,  yielding 


r  ""I1      X  1         .2 

f      +   /   ff  dn  +  3  /   (1  -  f   )  dn 
L   Ji-1    i-1  i-1 


'V.  (f  3XnT  "  f   8lnT)dn  "  ^2   '   at  dn 
l-l  *      a*   l-l 


i  9£nu 

^2   /  (f '  "  1)  —^  dn  =  0  (43) 

a*   i=l 


Expanding  equation  (4  3)  in  Taylor  series  and  applying 
equation  (36)  to  the  time  derivative  terms: 
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)    af '      _  faf" 


i-1 


=    tTofk   +    Vk-1    +   Vk-2^  <44> 


i  9£nu  9&nu 

f      (f '-   1)    — kt-S  dn   ■  — ST^  {l  f  I 
i-1  8t  8t         L     Ji-1 


6n>  [T   ku     +   T.ilnu  +T    Jlnu  ](45) 

°        ek        L        ek-l      Z        ek-2 


the    integrated  boundary    layer  equation   becomes: 


[ 


♦  f     (l+d0)f  +  (d1f1_1+d2fi.2)f  -^    UTofk+Tlfk.1+T2fk.2] 


-f[T  tou   +  T 
°    Sk 


1    ek-l    2    ek-2   Ji-1 
-(1  +  3  +  2dQ)  [f !!  XMX  +  f  V  XM2  +  fi"xM3  +  f;["1XM4] 


-2[fi  ZMX  +  f\  ZM2  +  f.  ZM3  +  fi_]_ZM4] 


■6n{6  +  -—7  [T  £nu   +  T,  £nu     +  T,£nu     ]  }  =  0   (46) 
a*2    °    ek     X    ek-l     ^    ek-2 


The  Taylor  series  expansions,  XM  and  ZM,  of  the  integrals 
in  equation  (43)  are  listed  in  Appendix  A. 

G.   RECURRENCE  FORMULAS 

The  boundary  layer  equations,  equation  (46),  and  the 
Taylor  series  expansions  of  the  primary  variables,  equations 
(27)  through  (29),  together  with  the  boundary  conditions, 
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equations  (23)  through  (26),  are  solved  by  Newton-Raphson 
iteration  of  4N  simultaneous  equations  in  4N  primary 
variable  unknowns.   These  are  summarized   as  follows: 


Primary  Variables 
f 

K 

f- 

No.  of  Eqns.  Req : 


No.  of  Equations 


4N 


Taylor  Series  Expansions 
Boundary  Layer  Equations 
Boundary  Conditions 
No.  of  Eqns.  Available: 


Equation  Nos 

(27)-(29) 

(46) 
(34)-(26) 


No.  of  Equations 

3(N-1) 
N-l 
4 


4N 


The  Newton-Raphson  solution  technique  is  illustrated  by 
considering  two  simultaneous  nonlinear  equations: 


F(x,y)  =  0 


G(x,y)  =  0 


(47) 


which  have  the  solution  x=X  and  y=Y.   Defining  x   and  y   as 

m     u  m 

the  values  of  x  and  y  for  the  m   iteration,  the  desired 
solution  f(X,Y)  is  expressed  in  the  Taylor  series  expansions 


0  =  F(X,Y)  =  F(x  ,y J  +  (X  -  x J 
mm  m 


3x 


(Y  "  V 


9y 


0  =  G(X,Y)  =  G(x  ,y J  +  (X  -  x J 
mm  m 


9x 


(Y  -  y  J 
m 


9G(x  ,y J 

mm 

9y 


(48) 
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In  the  Newton- Raphs on  method  (X,Y)  is  replaced  by 

^xm+l'ym+l^  and  the  nonlinear  terms  in  x  and  y  appearing  in 
equations  (4  8)  are  ignored  yielding  the  set  of  recurrence 
formulas : 

Axm  5S +  A^m  5y =  -F(xm'V 

8G(x  ,y  )        3G(x  .y) 

***  -t^  +  ^  —^-^  =  -G<=vV  <«> 

where : 

Axw  =  x^,  -  x  Ay^  =  y  .,  -  ym  (50) 

m    m+ 1    m  u  m   ■*  m+ 1   J  m 

The  Ax   and  Ay   are  corrections  to  be  added  to  x   and  y 

m      Jra  m     J  m 

respectively  yielding  x  .  -  and  y  ,,.   F(x  ,y  )  and  G(x  ,y  ) 

c  J  m+ 1     u  m+ 1      m  ■*  m         m  J  m 

are  the  values  of  the  original  functions  F(x,y)  and  G(x,y) 
evaluated  for  the  value  of  the  variables  during  the  m 

iteration.   Noting  that  F  (x  ,y  )  and  G(x  ,y  )  approach  zero 

mm         m  u  m   cir 

as  the  corrections  approach  zero,  it  is  appropriate  to  con- 
sider them  as  errors  associated  with  the  original  equations 

(47). 

Differentiating  the  Taylor  series  expansions,  equations 

(27)  through  (29)  ,  with  respect  to  the  primary  variables  in 

accordance  with  the  derivation  of  equation  (49)  yields  the 

following  recurrence  formulas  for  the  m    iteration: 

2         3         3 
(-1)  Afi+1+  (1)  Af±+  (6n)  tf[+  (^-)  Af  %  (£fl-)  Af  JV(^L-)  Af^—  ERROR 

(51) 
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2  2 

(-1)  Af  '   +(l)Af  !+(6n)  Af  "+(^-)Af  !"+(^-)  Af  ."'  =  -  ERROR   (52) 
1  + x       l        l    3     i     6     1+.L 


(-l,)Af^+1+(l)Af™  +  (^r).Af™  +  (^4*^+1  ="  ERROR        (53) 


where  Af.  represents  the  correction  for  f.,  the  coefficients 
of  the  corrections  represent  the  partial  derivatives  of  the 
Taylor  series  expansions  with  respect  to  the  primary  var- 
iables, and  the  ERRORS  are  equations  (29)  through  (31) 
respectively,  evaluated  for  the  values  of  the  variables 
during  the  iteration. 

The  recurrence  formulas  for  the  boundary  layer  equations 
are: 


Af   +  Af  (l+dQ)f  +  f(l+dQ)Af  +  Af  (d]_f£_1+  d2f£-2} 

— Hr  {T  Af  -  Af(T  Ami   +  T,£nu     +  T«£nu     )}i 
a*2    °  °    ek    1    ek-l    2    ek-2   Ji-1 


-2(1+  3  +  2dQ)  {AfiXM1+AfiXM2  +  Af±  XM3+  Af^^XM  } 


•2{AfiZM1  +  AfiZM2  +  Af.  ZM3+  Afi_1ZM4>=  -ERROR      (54) 


where  ERROR  is  the  left-hand-side  of  equation  (46)  evaluated 
for  the  values  of  the  variables  during  the  iteration. 

Finally,  the  recurrence  formulas  for  the  boundary  con- 
ditions, equations  (23)  through  (26),  are: 


Af '  =  -  ERROR  =  -  (f ' )  (55) 

w  w  m 


Af   =  -  ERROR  =  -  (f  )  (56) 

w  w  m  v   ' 
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Af   =  -  ERROR  =  1  -  (f  )  (57) 


Afe  =  -  ERROR  =  -  (fe)m  (58) 


Equation  (54)  may  be  further  refined  by  collecting  terms 
and  writing  the  equation  in  terms  of  correction  coefficients 
as  follows: 


ClAfi  +  C2Afi_1  +  C3Afi  +  C4Afi_1 


+  C5Afi  +  C6Af._,  +  C7Af.   +  C8Af . ,-  -  ERROR   (59) 


The  correction  coefficients  CI  through  C8  are  given  in 
Appendix  B. 

H.   MATRIX  FORMULATION 

The  coefficients  for  the  recurrence  formulas,  equations 
(51)  through  (53)  and  (55)  through  (59) ,  evaluated  for  the 
m   iteration  form  a  square  matrix  I  A J  with  4N  rows  and 
columns.   The  matrix  equation  is: 

A   I  AV  I  =  -   E  (60) 

L  JL  J   L  J 

where  the  column  matrix  I  AV J  is  the  matrix  of  corrections  to 
the  primary  variables,  and  the  column  matrix  E   is  the 
matrix  of  errors.   The  matrix  A  contains  many  zeros,  and 
great  savings  in  computer  storage  space  and  computation  time 
may  be  realized  through  proper  ordering  of  the  equations  and 
variables. 
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The  equations  are  first  divided  into  linear  and  nonlinear 
sets,  and  at  the  same  time,  the  variables  are  classified  as 
linear  or  nonlinear,  resulting  in  the  partitioned  matrix 
equation : 


(61) 


AL   !   BL 

AVL 

EL 

ANL  ;  BNL 

AVNL 

ENL 

The  division  of  variables  into  linear  and  nonlinear  sets  is 
generally  arbitrary  except  that  the  division  must  be  made  to 
insure  that  the  square  matrix  AL   is  not  singular.   Adopt- 
ing the  variable  groupings  of  Bartlett  and  Kendall  (Ref.  16) 
the  linear  and  nonlinear  corrections  are  arranged  as  follows: 


AVL(Af2  ,Af  3,  .  .  .  ,Afe,Afw,Af2,  .  .  .  ,  Af  e  ,  Afw,  Af  2  ,  .  .  .  /Af^) 
AVNL(Afw,Af^,Afj;Af2",...,Af") 


The  linear  and  nonlinear  equations  are  sequenced  within  the 
partitioned  square  matrix  :  A  j  as  shown  in  Figures  1  and  2, 
which  demonstrate  the  case  of  4  nodes  across  the  boundary 
layer. 

Expanding  equation   (61) ,  the  resulting  linear  and  non- 
linear equations  are,  respectively: 


F   IT    1 

AL    AVL  j 

j  ANL   |  AVL  j 


BL 


BNL 


AVNL 


AVNL 


I  =   -ENL 


(62) 
(63) 


Solving  equation  (62)  for  the  linear  corrections: 
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r  n_1~   i r      i    r   n~1  r     i 

AVL   =  -!  AL  !    BL  ;  ,  AVNL   +   AL  |    \    -EL 

L   j    L     J 


(64) 


r   -» 
Examination  of  the  '  AL  |  and  ,  BL  i  matrices  reveals  that 


IAT.1 

r   t  -r  —\ 

the  product  !  AL  ■      !  BL  |  may  be  performed  given  only  the  nodal 

spacing  and  must  be  determined  only  once  for  a  given  problem, 
thus  : 


AVL 


BA1  j  ;  AVNL 


+r 


ELA 


(65) 


where : 


BA1  I  =  i  AL 


TTblI 

j   L    J 


ELA 


1  -EL 


(66) 
(67) 


.T  i 

Applying  equation  (65)  to  equation  (63) ,  and  solving  for 
the  nonlinear  corrections: 


1       r 

AVNL    ;  =  I 


where : 


fH 


II 


BNL    • =       BNL    • -   j    ANL      ,    BAl    \ 

.j     L       J     L       Jl       j 


ENL 


=  i    -ENL 


i  r    ir    i 

j  -  |   ANL      |   ELA 


(68) 

(69) 
(70) 


I.   SOLUTION  PROCEDURE 

The  matrix  solution  proceeds  according  to  the  following 
algorithm: 

r    i"1 

1.  Given  the  nodal  spacing  determine  |  AL     and  form 
the  product  ["BAl") . 

2.  Given  initial  values  of  the  primary  variables  deter- 
mine r-EL~]  and  f~-ENLl,  and  form  the  product  j~ELA  |. 
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r -]      r -i  p    1 

3.  Form  !  BNL  ,  and   ENLj  from  the  coefficients  of  j  ANL  i 
and  j~BNL  ;'. 

4.  Invert  j  BNL  ;  and  form  the  product   BNL j    ENL  ', 
which  is  the  solution  for  j  AVNL~j. 

5.  Solve  for  [avlJ  from  equation  (64)  using  I  AVNLJ 
as  determined  i4i  step  4 . 

6 .  Add  the  linear  j  AvJ  to  the  corresponding  primary 
variables  to  complete  the  iteration. 

r     i       r      i 

7.  Compute  the  new  errors  . -EL J  and  j -ENL  j,  and 

test  for  a  maximum  allowable  error.   !£f  thi's  criterion 
has  not  been  reached,  return  to  step  2  using  these 
errors  to  commence  the  new  iteration. 
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III.   RESULTS  AND  DISCUSSION 

The  computer  program  has  been  used  to  investigate  three 
basic  categories  of  boundary  layer  problems  using  a  flat 
plate  model;  steady  Blasius  flow  (zero  pressure  gradient) , 
steady  Howarth  flow  (uniform  adverse  pressure  gradient) , 
and  oscillating  Blasius  flow.   The  velocity  profiles  which 
resulted  from  these  solutions  are  compared  to  previously 
reported  results  and  presented  in  Figures  3,  4,  and  5. 

A.  STEADY  BLASIUS  FLOW 

In  Figure  3,  computed  values  corresponding  to  three 
nodal  spacing  choices  are  superimposed  on  the  Blasius  pro- 
file.  Agreement  between  the  numerical  solution  and  the 
classical  solution  is  better  than  0.08  percent  for  as  few 

as  10  nodal  points. 

. 

B.  STEADY  HOWARTH  FLOW 

Results  for  a  steady  flow  in  a  uniform  adverse  pressure 
gradient  are  compared  with  the  solution  of  Howarth  (Ref .  11) 
in  Figure  4.   The  comparison  of  the  computer  solution  with 
that  of  Howarth  is  accomplished  by  correlating  the  dimen- 
sionless  parameter  x*  (=b,s/U  ),  which  characterizes 
Howarth * s  velocity  profiles.   Agreement  between  the  two 
solutions  is  generally  good  with  differences  no  greater  than 
4.35  percent  for  22  nodal  points.   The  differences  encoun- 
tered may  be  reduced  somewhat  by  increasing  the  free  stream 
velocity  of  the  test  problem.   This  results  in  larger 
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FIGURE  3.  STEADY  BLASIUS   FLOW  DATA  COMPARISON 
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stream-wise  pressure  differences  in  the  equations  of  motion 
and  consequently  lesser  roundoff  errors. 

C.  OSCILLATING  BLASIUS  FLOW 

The  unsteady  capabilities  of  the  numerical  procedure 
were  investigated  in  the  solution  of  Blasius  flow  with 
sinusoidal  oscillations  superimposed  on  the  mean  flow. 
Utilizing  the  critical  oscillation  frequency  go   (=.6U  /s, 
for  a  flat  plate)  of  Lighthill's  analysis  (Ref.  2)  to 
determine  the  solution  approximation  with  which  to  compare 
the  numerical  results,  the  profile  of  velocity  fluctuations 
which  resulted  was  compared  with  the  "quasi-steady" 
analysis  as  shown  in  Figure  5.   Agreement  between  the  calcu- 
lated velocities  and  the  Lighthill  profile  is  within  15 
percent  deep  in  the  boundary  layer  and  within  5.5  percent 
in  the  outer  regions.   Noting  that  Lighthill's  analysis  is 
based  on  linearized  small  perturbation  theory  and  that  the 
amplitude  of  oscillation  in  the  test  problem  of  0.1  U   can- 
not be  considered  a  small  perturbation,  better  agreement 
could  scarcely  be  expected. 

D.  NUMERICAL  PROCEDURE 

The  computer  program  has  been  operated  on  both  IBM  360/67 
and  CDC  6600  computer  installations.   Core  storage  require- 
ments are  130,000  bytes  for  a  maximum  array  size  encompassing 
25  nodal  points  across  the  boundary  layer,  50  streamwise 
stations  and  50  time  increments  per  oscillation  cycle.   A 
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typical  25  nodal  point  Blasius  solution  was  performed  in 
7  seconds  on  the  CDC  6600. 

Limitations  in  the  present  numerical  capabilities  of 
the  program  precluded  a  more  extensive  investigation  of 
parametric  effects,  however  the  results  qualitatively 
demonstrated  the  feasibility  of  the  computational  procedure 
Further  refinements  to  the  computer  program  should  enable 
a  more  extensive  investigation  to  include  quantitative  com- 
parisons with  experimental  results. 
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IV.   CONCLUSIONS  AND  RECOMMENDATIONS 

The  feasibility  of  the  Unsteady  Boundary  Layer  Integral 
Matrix  Procedure  in  the  treatment  of  simple  steady  and 
unsteady  flow  problems  has  been  demonstrated.   Favorable 
comparison  with  classical  results  provides  a  firm  base  for 
continued  investigation  in  the  unsteady  regime. 

Continued  research  utilizing  the  program  should  include 
the  following  revisions: 

1.  Alteration  of  the  present  time  derivative  component 
arrays  (FTM1)  and  FTM2)  to  insure  that  all  values  of  the 
components  are  stored  for  each  nodal  point,  instead  of 
retaining  only  the  values  at  the  wall. 

2.  Incorporation  of  Blasius  profiles  as  first  estimates 
for  the  primary  variables  in  order  to  decrease  the 

time  required  to  complete  the  first  nodal  iteration 
procedure. 

3.  A  study  of  array  usage  aimed  at  reducing  the  re- 
quired storage  and  subsequent  conversion  to  double 
precision  computation  in  order  to  minimize  round-off 
interference. 

4.  Addition  of  the  capability  of  handling  compressible 
flows  in  order  that  higher  velocities  may  be  introduced. 

5.  Inclusion  of  turbulent  boundary  layer  capabilities 
perhaps  through  the  use  of  eddy  transport  properties. 
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APPENDIX  A 
TAYLOR  SERIES  EXPANSIONS  OF  INTEGRALS  (Ref.  16) 
i  r    Hi      i 

/     ff"  dn  =  i  ff '  j       -  /     f       dn  (Ai) 

i-1  L        Ji-1      i-1 


/       f'2    dn    =    f;J_XM1    +    f^XM2    +    f!UXM3    +    fl'^XM^  (A2) 

xmx  .  «n  (    f :  -  fj  .^L  +  f »'  i£  ,  f-:i  *£,       (A3) 

XM2=W(|f:-f^      .fi-igj+f^^  (A4) 

xm,  =  6n3£  f!   -  f'.'-^  +  f!"J3*£  +  f!"     Sfinlj      (A5) 

3  v8      i  i  120  i     420  1_1     504 

xm,  =  6n3(^-f:    -  f?  il    +  f-W       f.«      6ni} 

4  v24    i  i    30  i      504  1-1    252; 
^^TOf  ^"aM^    =    -[doff'    +   dlVlf'    +    d2f£-2f,]i-l 


+2dQ[f^XM1    +    f^XM2    +    fi"XM3    +    fi-lXM4] 

+2[f!zJYL     +    f'.'zM-    +    f!"  ZM0    +    f!" ,]  (A7) 
l       1            l       2            l         3            l-l 

2  2 

ZM±    =    6n     (       YMX    -    YM2    iJ2-  +    YM3    ^-      +    YM4  ^_)                      (A8) 

2  2 

ZM2    =-6n2(|  YMX    -    YM2    ^-  +    YM3    ijjgj-  +    YM4  ^-)                      (A9) 

ZM3    -    6n3(i  YM,    -    YM,^  +    YM      W    +    ^  5^,                (A1Q) 
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ZM3=    6n3(^YM1-YM2^+YM3^+YM4|^ 


(All) 


YM,    =    d,f0    ,     •    +    d0f„    0     . 
1  1    £-l,i  2    £-2,1 


(A12) 


YM„    =    d,f„     ,     .     +    d„f0    _     . 
2  1    £-1,1  2    £-2,i 


(13) 


YM3   =    dlfi-l,i    +    d2f£-2,i 


YM4    =    dlf£-l,i-l    +   d2f£-2,i-l 


(A15) 
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APPENDIX  B 
BOUNDARY  LAYER  EQUATION  CORRECTION  COEFFICIENTS 


CI  =  (1+d  )f!-£±L.  {T  (1-inu    )-(T,£nu     +  T_£nu     )}   (Bl) 
o   1    j.z    o       e     1    e,  ,    2    e,  ~ 
a*  k-1        k-2 


C2  =-(l+d  )f.'  +  -£-H_  {T  (l-£nu  )-(Tn£nu     +  TJnu     )}   (B2) 
°  1"1  a*2   °  1    ek-l    2    ek-2 


C3  =  (l+do)fi+(dlfjl.    +d2f£_    )-2(l+6+2do)XM1-2ZM1     (B3) 


C4  =-(1^)^  "  (^Vi,!-!  +  V£-2,i-l>  (B4) 


C5  =  1  -  2 (l+3+2dQ)XM2  -  2ZM2  (B5) 


C6  =  -1  (B6) 


C7  =  -2(l+3+2d  )XM   -  2ZM3  (B7) 


C8  =  -2(l+3+2d  )XM„  -  2ZM.  (B8) 

o    4       4 


43 


APPENDIX  C 


PROGRAM  DESCRIPTION 


A.   GENERAL 

The  computer  program  referred  to  as  UBLIMP,  is  written 
in  FORTRAN  IV  source  language  and  has  been  operated  on  both 
IBM  360/67  and  CDC  6600  computer  installations.   Core 
storage  requirements  are  130,000  bytes  for  a  maximum  array 
size  encompassing  25  nodal  points  across  the  boundary  layer, 
50  streamwise  stations  and  50  time  increments  per  oscillation 
cycle . 

The  program  consists  of  a  main  program  (UBLIMP)  and  21 
subroutines  which  are  divided  into  two  sections  as  shown 
in  Figure  6.   LINK  0,  common  to  both  sections,  contains  the 
main  program  and  certain  service  routines  which  are  common 
to  both  sections  in  the  fully  expanded  program.   LINK  1 
sets  up  the  boundary  conditions  and  LINK  2  iterates  to  solve 
the  boundary  layer  equations. 

The  program  may  be  applied  to  steady  or  unsteady  flow, 
two-dimensional  cartesian  or  axisymmetric  bodies  with  similar 
or  nonsimilar  profiles.   It  allows  quadratic  or  cubic  curve 
fits  of  the  primary  variables.   Program  options  are  con- 
trolled by  the  control  variable  KR  and  are  detailed  in  the 
input  instructions  for  the  UBLIMP  program  in  Appendix  D. 
A  program  linkage  schematic  is  shown  in  Figure  7  and  a 
functional  flow  chart  is  shown  in  Figure  8.   Correlation 
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LINK  1 


•  UBLIMP 
—  ZEROIT 


—  MATONE 


LINK  0 

■—  SLOPQ 

—  TAYLOR 

—  RERAY 

—  SETUP 

— ITERAT 

—  RECASE 

—  LINCER 

—  INPUT 

—  ABMAX 

—  STDETA 

—  NONCER 

—  LINMAT 

LINK  2 

—  IMONE 

—  FIRSTG 

—  AMSET 

■—  REFCON 

—  I ONLY 

—  HISTXI 

—  OUTPUT 

FIGURE  6.   COMPUTER  PROGRAM  OVERLAY  STRUCTURE 
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CASE 


A    READ  INPUT  DATA  FOR  NEW  CASE  WHICH  DIFFERS  FROM  THAT 
OF  PREVIOUS  CASE. 


B     IF   NODAL  SPACING  DIFFERENT  FOR  THIS  CASE  SET  UP  MATRICES 
FOR  TAYLOR  SERIES   EXPANSIONS  AND  LINEAR  BOUNDARY  CON- 
DITIONS  FROM  NODAL  POINT  DISTRIBUTION.      INVERT   MATRIX 
OF  COEFFICIENTS   OF  LINEAR   CORRECTIONS     (THIS,  IN   EFFECT, 
EXPRESSES    LINEAR  CORRECTIONS  IN   TERMS  OF   NONLINEAR 
CORRECTIONS). 


C    COMPUTE  OR  READ  FIRST  GUESSES  FOR  PRIMARY  VARIABLES 
OR  USE   PREVIOUS    VALUES. 


STATION 


D    CALCULATE   BOUNDARY    LAVER    EDGE   CONDITIONS   AROUND 
THE   BODY. 


TIME 


E    COMPUTE   NONSIMILAR    TERMS     (THOSE  INVOLVING    STREAM- 
WrSE    DERIVATIVES). 


ITERATION 


F    EVALUATE    LINEAR   ERRORS. 


SPLINE 
POINT 


G    EVALUATE    ERRORS   AND  COEFFICIENTS  FOR   NONLINEAR 
EQUATIONS    AND  CORRECT    NONLINEAR     ERRORS    FOR 
LINEAR     ERRORS. 


JH    INVERT    MATRIX   OF  NONLINEAR   CORRECTION  COEFFICIENTS. 

I 


I     EVALUATE    NONLINEAR    AND    LINEAR  CORRECTIONS. 
CORRECT     PRIMARY  VARIABLES   AND   TEST  ERRORS   FOR 
CONVERGENCE. 


E 


PKIN-    OLTPUT    DATA  FOR  CURRENT    STATION    AND  TIME. 


K    INITIALIZE    TIME. 


STOP 


FIGURE  8.   FUNCTIONAL    FLOW  CHART  OF  UBLIMP  PROGRAM 
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between  the  two  is  indicated  by  letters  running  from  A 
through  K. 

B.   PROGRAM  FUNCTIONS 

1.  Main  Program  (UBLIMP) 

The  main  program  provides  linkage  between  the 
boundary  conditions  section  and  iterative  solution  section 
of  the  program.   Initial  core  zeroing  and  the  call  for 
solution  output  are  also  initiated  by  the  routine. 

2.  Subroutine  ZEROIT 

ZEROIT  performs  the  initial  zeroing  of  all  named 
common  blocks  to  eliminate  the  possibility  of  introducing 
random  numbers  in  the  solution. 

3.  Subroutine  SETUP 

SETUP  is  the  control  routine  for  setting  up  the 
boundary  layer  edge  conditions  and  the  streamwise  derivatives 
for  a  new  time,  station  or  case. 

4 .  Subroutine  RECASE 

RECASE  is  the  control  routine  for  the  input  of 
boundary  layer  data  and  initializes  control  variables  for 
the  treatment  of  surface  discontinuities. 

5.  Subroutine  INPUT 

INPUT  reads  all  the  data  for  a  problem  into  the 
computer. 

6.  Subroutine  STDETA 

STDETA  permits  simplification  of  the  input  data  deck 
by  providing  standard  nodal  point  distributions  should  the 
user  not  care  to  provide  one. 
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7.  Subroutine  LINMAT 

LINMAT  sets  up  the  matrices  for  the  Taylor  series 
expansions  and  linear  boundary  conditions  from  the  eta 
spacing  and  solves  to  express  the  linear  corrections  in 
terms  of  the  nonlinear  corrections. 

8.  Subroutine  MATONE 

MATONE  performs  operations  on  a  column  of  a  matrix 
of  coefficients   BL   or  on  a  column  of  a  matrix  of  errors 

-EL   (designated  E  in  the  subroutine  call  list)  such  as  to 
form  TalI  [~E  1. 

9.  Subroutine  FIRSTG 

FIRSTG  computes  first  estimates  of  the  primary 
variables  based  on  the  nodal  point  distribution  should 
the  user  not  care  to  provide  them. 

10.  Subroutine  REFCON 

REFCON  calculates  the  boundary  layer  edge  conditions 
and  sets  up  the  wall  boundary  conditions. 

11.  Subroutine  SLOPQ 

SLOPQ  defines  a  cubic  equation  for  a  set  of  points, 
calculates  the  average  slope  at  each  point  and  integrates 
the  cubic  equation  between  each  pair  of  points. 

12.  Subroutine  HISTXI 

HISTXI  computes  terms  involving  the  axial  deriva- 
tives and  time  derivatives  and  stores  those  upstream  quan- 
tities needed  for  these  difference  relations. 
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13.  Subroutine  TAYLOR 

TAYLOR  calculates  the  coefficients  in  the  Taylor 
series  expansions  of  the  integrals  which  appear  in  the 
integral  form  of  the  boundary  layer  equations. 

14.  Subroutine  ITERAT 

ITERAT  is  the  control  program  for  performing  the 
boundary  layer  iteration. 

15.  Subroutine  LINCER 

LINCER  evaluates  errors  for  the  linear  equations 
and  determines  the  maximum  linear  error. 

16.  Subroutine  ABMAX 

ABMAX  searches  a  given  array  for  the  entry  with 
the  maximum  absolute  value. 

17.  Subroutine  NONCER 

NONCER  performs  that  portion  of  a  boundary  layer 
iteration  having  to  do  with  the  solution  of  the  nonlinear 
corrections  of  the  nonlinear  equations,  computes  the  damp- 
ing factor  and  applies  it  to  the  corrections  and  then 
applies  the  damped  corrections  to  the  primary  variables. 

18.  Subroutine  IMONE 

IMONE  evaluates  the  coefficients  of  (i-1)  cor- 
rections for  the  i    nonlinear  equations,  where  i  is  the 
i    nodal  point  in  the  boundary  layer. 

19.  Subroutine  AMSET 

AMSET  calculates  the  contributions  of  the  I  ANlJ  and 


applies  them  to  j  BNL   and   ENL  . 
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20.  Subroutine  IONLY 

IONLY  evaluates  the  coefficients  of  the  i  cor- 
rections for  the  i  nonlinear  equations,  where  i  is  the 
i    nodal  point  in  the  boundary  layer. 

21.  Subroutine  RERAY 

r  "I 

RERAY  replaces  a  matrix  I  C  j  with  its  inverse  and 

r  "i 
forms  the  product  of  the  inverse  with  another  matrix  D  ,. 

22.  Subroutine  OUTPUT 

OUTPUT  prints  a  standard  boundary  layer  output  block 
for  a  converged  solution,  or  at  the  end  of  each  iteration 
if  desired. 
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APPENDIX  D 
INPUT  INSTRUCTIONS  FOR  THE  UBLIMP  PROGRAM 


GROUP  1   CONTROL  CARD,  TITLE 

CARD  1   FORMAT (20 II, 15A4) 

FIELD  1   (Columns  1-20)   This  is  the  variable  KR  which  is 
used  to  control  the  various  program  options. 

COLUMN  1   Determines  whether  a  new  set  of  ETA  values 
is  to  be  input  for  the  present  case. 

0  Uses  resident  values  from  previous  case. 

1  Values  input  by  user  (mandatory  for  first  case) . 

COLUMN  2   Designates  the  type  of  first  guesses  to  be 
utilized  for  primary  variables. 

0  Uses  built-in  relations  to  calculate  first 
guesses . 

1  First  guesses  input  by  user. 

2  Uses  resident  values  from  previous  case  (cannot 
be  used  for  first  case) . 

COLUMN  3   Determines  the  treatment  of  streamwise 
derivatives . 

0  Performs  similar  solution  at  each  streamwise 
station. 

1  Considers  two-point  difference  relations  at  all 
stations  except  that  a  similar  solution  is  per- 
formed at  the  first  station  for  non-blunt  bodies 
and  at  the  first  two  stations  for  blunt  bodies. 

2  Considers  three-point  difference  relations  at  all 
stations  except  that  a  similar  solution  is  per- 
formed at  the  first  station  and  a  two-point  solu- 
tion is  performed  at  the  second  station  for 
non-blunt  bodies;  similar  solutions  are  performed 
at  the  first  and  second  stations  and  a  two-point 
solution  is  performed  at  the  third  station  for 
blunt  bodies;  and  a  two-point  solution  is  per- 
formed for  the  first  station  after  a  discontinuity. 
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COLUMN  4   Determines  when  output  block  is  to  be 
printed. 

0  Output  block  printed  for  converged  solution  or 
for  nonconverged  solution  after  50  iterations. 

1  Output  block  printed  after  each  iteration. 

COLUMN  5  Determines  treatment  of  Entropy  Layer. 
(Not  used  in  present  program.  KR(5)=4 
is  input  to  maintain  input  data  sequencing.) 

COLUMN  6   Designates  body  shape. 

0  Axisymmetric  blunt  body. 

1  Planar  blunt  body 

2  Axisymmetric  sharp  body. 

3  Planar  sharp  body. 

4  Axisymmetric  or  planar  shape  which  has  no  sharp 
tip  or  blunt  stagnation  point,  such  as  a  nozzle. 

COLUMN  7   Steady,  zero  pressure-gradient  noise  filter. 

0  Normal  computations  of  DUDS. 

1  DUDS  is  forced  to  zero. 

COLUMN  8  Designates  form  in  which  wall  mass  fluxes 
are  input.  KR(8)  is  not  utilized  if  wall 
mass  fluxes  are  not  input. 

0  Wall  fluxes  input  in  LBS/SEC  FT**2 

1  Wall  fluxes  input  in  normalized  form  (divided  by 
-ALPHASTAR) . 

COLUMN  9   Together  with  KR(ll)  this  designates  the 
type  of  wall  boundary  conditions. 

0  Assigned  stream  function  at  the  wall. 

1  Assigned  mass  flux  at  the  wall.   KR(9)=1  is  used 
if  zero  mass  flux  is  used. 

COLUMN  10  Determines  the  type  of  curve  fits  employed 
to  represent  the  primary  variables. 
(KR(10)=0  is  recommended  for  accuracy  for 
most  problems,  however  for  severe  problems 
KR(10)=1  is  better  since  the  cubics  can 
become  poorly  behaved.) 
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0  Utilizes  connected  cubics. 

1  Utilizes  connected  quadratics  except  for  the 
outermost  segment  where  connected  cubics  are 
used. 

COLUMN  11   Together  with  KR(9)  this  designates  the 

type  of  wall  boundary  conditions.  (Presently 
KR(11)=0  is  the  only  active  value  and 
requires  an  assignment  of  wall  temperature.) 

COLUMN  12   Provides  access  to  standard  nodal-point 
distributions  built  into  the  program. 

0  User  inputs  data  for  nodal-points. 

1  8  points  with  more  concentration  at  the  wall. 

2  8  points  equally  spaced. 

3  10  points  equally  spaced. 

4  12  points  equally  spaced. 

5  14  points  equally  spaced. 

6  15  points  with  more  concentration  at  the  wall. 

7  18  points  equally  spaced. 

8  22  points  equally  spaced. 

9  25  points  equally  spaced. 

COLUMN  13   Permits  the  assignment  of  a  convergence 
damping  factor.   (This  factor  is  over- 
ridden if  a  smaller  damping  factor  is 
computed  internally  by  some  constraint.) 

0   No  damping  factor  is  assigned. 

J   If  J  is  greater  than  zero,  corrections  are 
damped  uniformly  by  J/10. 

COLUMN  14   Non-zero  entry  causes  a  complete  set  of 

primary  variables  to  be  output  for  future 
use  as  first  guess  inputs. 

COLUMN  15   Non-zero  entry  provides  debug  output  for 
first  guesses  and  linear  matrices. 
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0  No  debug. 

1  First  guesses  are  dumped. 

2  Linear  matrices  before  and  after  inversion  are 
also  dumped. 

COLUMN  16   Non-zero  input  indicates  the  pressure  pro- 
file has  been  input  in  the  form  of  pressure 
coefficients . 

COLUMN  17   Non-zero  entry  provides  debug  output  for 
coefficients . 

0   No  debug. 

Y   For  (Y+l-ITS)  greater  than  zero,  where  ITS  is  the 
number  of  the  current  boundary  layer  iteration, 
the  coefficients  which  combine  to  make  up  the 
nonlinear  equations  are  dumped  and  the  derivatives 
of  the  nonlinear  equations  with  respect  to  the 
nonlinear  variables  are  dumped  before  and  after 
inversion . 

COLUMN  18   Not  used. 


COLUMN  19 


Non-zero  entry  provides  debug  output  for 
linear  and  nonlinear  errors. 


COLUMN  2  0   Determines  which  of  the  two  internal  sets 

of  first  guesses  are  to  be  used  if  the  user 
does  not  provide  them. 

0  Straight  line  Blasius  inputs  are  calculated  for 
the  stream  function  and  velocity  ratio,  and  zero 
is  input  for  the  first  and  second  derivatives  of 
the  velocity  ratio. 

1  Straight  line  Blasius  inputs  are  calculated  for 
all  the  primary  variables. 

FIELD  2   (Columns  21-80)  CASE.   Title  of  the  case 

(ALPHA-NUMERIC) .   Used  for  identification  of 
printed  output. 

GROUP  2   TIMES  AND  STATIONS 

CARD  1   FORMAT (13) 

FIELD  1  (Columns  1-3,  right- justified)  NITEM  Number  of 
time  grid  points  per  cycle  when  considering  an 
unsteady  solution,  otherwise  input  1.  (max=50) 
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CARD  2   FORMAT (E10. 4) 

FIELD  1   (Columns  1-10)  TIME(l)   Value  to  be  used  in 

the  output  block  to  identify  a  time  of  solution. 

CARD  3   FORMAT (13) 

FIELD  1   (Columns  1-3,  right- justified)  NS.  Number  of 
streamwise  stations.   (max=50) 

CARD  4   FORMAT (8E10. 4) 

FIELD  1,2 ,...8  per  card (10  columns  per  field)  S(IS) 
Streamwise  distance  upon  which  the  boundary 
layer  solution  is  based  in  feet.   Blunt  body 
problems  should  start  with  S=0.0.   Sharp  body 
or  nozzle  problems  must  start  with  some  finite 
value  of  S.   The  boundary  layer  is  assumed 
similar  up  to  and  including  this  first  station. 
A  negative  entry  for  S  indicates  a  discon- 
tinuity at  that  station,  and  alters  the  differ- 
encing scheme  for  axial  derivatives. 

GROUP  3   NODAL  DATA   (Skip  this  GROUP  for  KR(12)^  0) 

CARD  1   FORMAT (13) 

FIELD  1   (Columns  1-3,  right- justified)  NETA  Number  of 
nodal  points  across  the  boundary  layer, 
(maximum  of  25) 

CARD  2   FORMAT (8E1 0.4) 

FIELD  1,2,... 8  per  card (10  columns  per  field)  ETA(I) 

ETA  stations  across  the  boundary  layer  starting 
at  the  wall  (ETA=0.0).   It  is  recommended  that 
a  value  between  5.0  and  6.0  be  assigned  at  the 
boundary  layer  edge.   UBLIMP  values  of  ETA  are 
equivalent  to  Blasius  values  divided  by  the 
square  root  of  two. 

GROUP  4   BODY  SHAPE  DATA 

CARD  1   FORMAT (2E1 0.4)  USED  ONLY  IF  KR(6)=0  or  1. 

FIELD  1   (Columns  1-10)  CONE.   Cone  half  angle  in  sphere- 
cone  shape  bodies.   Leave  blank  for  other  body 
shapes . 

FIELD  2   (Columns  11-20)  RNOSE.   Effective  nose  radius 
in  feet.   Used  to  calculate  stagnation  point 
velocity  gradient  from  Newtonian  relations. 
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CARD  2   FORMAT (8E10. 4)  USED  ONLY  IF  KR(6)=0,  2  or  4 

FIELD  1,2 ,...8  per  card (10  columns  per  field)  ROKAP(IS) 
This  is  the  local  body  radius  in  feet  normal 
to  the  body  centerline  raised  to  the  n  power, 
where  n  is  unity  for  axisymmetric  bodies  and 
zero  for  planar  bodies.   Therefore,  ROKAP  is 
unity  for  planar  bodies  and  the  local  body 
radius  for  axisymmetric  bodies.   For  planar 
bodies,  this  card  set  is  used  only  if  KR(6)=4. 
(A  special  input  format  can  be  used  for  spher- 
ical-nosed bodies  as  follows.  Set  ROKAP (1) 
equal  to  minus  the  nose  radius.   The  nose  radius 
is  then  set  to  -ROKAP (1)  and  ROKAP (1)  is  set  to 
zero.   If  subsequent  ROKAP  are  input  as  zeros, 
the  program  computes  ROKAP  from  S  for  a  spher- 
ical nose.   The  first  non-zero  entry  is  then 
ROKAP  at  the  tangent  point.   If  this  again 
followed  by  zeros,  linear  interpolation  is  used 
to  the  next  non-zero  entry  to  yield  ROKAP  along 
a  conical  afterbody.) 

GROUP  5   STAGNATION  DATA 

CARD  1   FORMAT (E10. 4) 

FIELD  1   (Columns  1-10)  PTET.  Local  stagnation  pressure 
in  atmospheres. 

GROUP  6   FIRST  GUESS  INFORMATION 

CARD  1   FORMAT (5 (2XE14.7) ) 

FIELD  1,2, ...5  per  card(16  columns  per  field)  F(I,J). 
First  guesses  for  stream  function  (F(1,J)), 
velocity  ratio  (F(2,J)),  shear  function 
(F(3,J)),  and  derivative  of  shear  function 
(F(4,J))  across  the  boundary  layer. 

GROUP  7   STREAMWISE  DISTRIBUTIONS  FOR  EDGE  CONDITIONS 

CARD  1   FORMAT (5 (2XE14.7) ) 

FIELD  1,2,... 5  per  card (16  columns  per  field)  PRE (IS). 
Ratio  of  local  static  to  stagnation  pressure 
for  the  mean  flow  condition.   When  KR(16)  is 
non-zero  PRE  is  input  in  the  form  of  negative 
pressure  coefficients  referred  to  PTET. 

CARD  2   FORMAT (Blank)  This  card  is  inserted  to  maintain 

data  continuity.   In  more  complex  formulations 
this  card  corresponds  to  a  control  flag  for 
updating  the  pressure  profile  during  the  com- 
putation . 

57 


GROUP  8   STREAMWISE  DISTRIBUTIONS  FOR  WALL  INPUT  CONDITIONS 

CARD  1   FORMAT (E10. 4) 

FIELD  1   (columns  1-10)  T.  Wall  temperature  in  degrees 
R  (for  the  present  this  is  the  free  stream 
static  temperature) . 

CARD  2   FORMAT (8E10. 4)  USED  ONLY  FOR  KR(9)=0 

FIELD  1,2,... 8  per  card(10  columns  per  field)  FW(IS). 
Wall  stream  function  (negative  for  mass 
addition) . 

CARD  3   FORMAT (Blank)   USED  ONLY  IF  CARD  2  IS  USED.   This 
card  is  inserted  to  maintain  data  continuity. 
In  more  complex  formulations  this  card  cor- 
responds to  a  control  flag  for  updating  the 
wall  stream  function  profile  during  the  compu- 
tation. 

CARD  4   FORMAT (8E1 0.4)   USED  ONLY  FOR  KR(9)=1 

FIELD  1,2 ,...8  per  card (10  columns  per  field)  RHOVW(IS) 
Total  mass  flux  at  the  wall  (LB/SEC  FT**2  or 
dimensionless  for  KR(8)=0  or  1  respectively). 
These  values  are  positive  for  mass  injection. 

CARD  5   FORMAT (Blank)   USED  ONLY  IF  CARD  4  IS  USED.   This 
card  is  inserted  to  maintain  data  continuity. 
In  more  complex  formulations  this  card  cor- 
responds to  a  control  flag  for  updating  the 
total  mass  flux  profile  during  the  computation 

GROUP  9   FREE  STREAM  DATA 


CARD  1   FORMAT (3E10. 4) 

FIELD  1   (Columns  1-10)   UMFS .  Mean  free  stream  velocity 
in  FT/SEC. 

FIELD  2   (Columns  11-20)  FREQ.  Frequency  of  mean  free 
stream  oscillations  in  HZ. 

FIELD  3   (Columns  21-30)  EPS.  Ratio  of  the  amplitude  of 
the  free  stream  oscillation  to  the  local  mean 
free  stream  velocity. 

LAST  CARD  FORMAT (Al)   The  purpose  of  this  entry  is  to  per- 
mit a  test  on  whether  or  not  a  new  case  is  to 
follow.   In  the  event  a  case  does  not  converge 
in  the  alloted  number  of  iterations,  any  re- 
maining cards  for  that  case  are  read  and  then 
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ignored  until  either  a  period  or  a  comma  is 
encountered  in  column  1.   A  comma  in  column  1 
indicates  that  another  case  is  to  follow,  while 
a  period  in  column  1  indicates  that  this  is  the 
last  card  in  the  input  deck. 

The  preceding  input  instructions  should  permit  the  user 
to  fully  exercise  the  UBLIMP  program  with  a  minimum  of 
difficulty.   For  more  detailed  information  on  various  rou- 
tines the  reader  is  referred  to  reference  18,  which  is  the 
user  manual  prepared  for  the  BLIMP  program. 
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COMPUTER  PROGRAM 

C  UNSTEADY    BOUNDARY    LAYER    INTFGPH    MATRIX    PROCEDURE 

f  (UPLTMD) 

REAL    Mt) 

COMMON/BUMCOM/BUMP,CORMA,EASE»ICORM,WDOTf I777, ISD, IX 

COMMON /FpGCOv/PEM( 50)  ,PTF(50)  ,DUES,UE<50,50)fPE(50,50) 
1,RHO,MU,T,DSIP(50)  ,IDSI  P 

CnMMON/INTCOM/KR(2  0)  ,K!N,K0'JT,MAT1I  .MAT  1  J.N FT  A.  I,  IS,  MS 
l,IT,NTIMF,NAM,NLEO,NNLFQ,NRNL,ITS,CASE(15  )  «  B(  8)  ,  MMEf  NO 
2N,*TFPf  NTTEi*,KP17  ,NBT,NBT2  ,  I  DENT  ,IJMFS  .  FRFO,  E^S 

rhMMON/poycnM/TiMF  (50)  »°RE(^0  )  , PTET,*  (  5  0  ) ,  ROK  AP  (  50  ) ,  PN 
10  SE  ,  VKAP.NDISCIPI  CC(50)  ,NS0(1Q  )  ,  MS  0(10)  ,  I  °RC  »R  AONO,  CO 
2NE,DELT, CYCLE 

COMMON/PR PCOM/PR.SC.XM (5) 

COMMON/VARCOM/F(5,25) ,FTMi (5 , 25 ) ,FTM2(5 , 25) 

CALL    7FPOIT 

1  FORWAT(Al) 

2  F0PMAKF1C.4) 

3  FQRMATf 1H1,///////////////5X10HRUN  TIME  =  FR.2,8H  SFCO 
1N0<  ) 

DATA  IAST/1H,/ 

D*t*  L*  FT  /1H.  / 

KIN=4 

Kni)T=P 

R(  1)  =  .5 

R( 2)=. 333333333 

R( 3)  =  . 166666666 

8(4)*. 125 

R( *)=.C41666666 

8(M  =  .C33333333 

B(7)=. 013888888 

8( 8)=. 0039*825* 

PHP=.C7651 

*U=3.71cE-07 

IT=1 

&e   mwf=-i 

EA*F=1. 

I?77=C 

!*=1 
*1  TTFM=l 
*2  CALL  SETU° 
*3  C^LL  ITFRfiT 

CALL  OUTPUT 

IF(N0N)43f44,4C 
to,     iTCM=ITfM+l 

no    4C    T=lfNETA 

00    4Q    J=l,5 

FTM2( J,I)=FTMl( J,I) 

FTM1( J,I)=F{ J, I ) 

4<5    CONTINUF 

IF(  TTFV-NITM)    42,42,45 

m  rs=is+i 

IF(IS-NS)    41,41,40 
40    RE^D(KIN?1)    JAST 

TF(  IA  ST-JAST)  47,46,^7 
47  TF(LAST-JAST)  ^0,48,40 
ca    «;TOP 

ENn 

SUBROUTINE    7FR0IT 
COMMON /BUMCOM/BUMCO(  8  ,1  ) 

COMMON /EDGCOM/EDGCO( 5 ,1)  ,EnGC(* ,10) ,EDG(1, 100), ED (5, 10 
100) 

COMMON /FPRCOV/FPPC 0(5,1)  ,ERRC (2  ,  10 ) , ERP ( 3, 100 ) 
COMMON/ETACOM/ETACO  <<=,!)  ,ETAC  (8  ,10  )  ,  ETA  (  2,  100  ) ,  ET  (  2,  10 
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ICO) 
COMMON /NT  ^CPM/HIECQ(6,1  )  ,HI<^C(3  ,10)  ,HIS  K,  100  ) 
COMMON /!MTCOV/INTCn(*  tl)  ,INTC(7  ,10) 
COMMON /NTNCOV/NONC (2,10) ,NON(8,100) 
COMMON /P«?  Mf  QM/pRMCQ(«,l)  ,PPMC  (7  ,10)  f  ORM(2,100  ) 
CnMMDN/PRPC0M/PRPCnf7,l) 

COMMON  /T^MCr3w/JFMCO(c  ,  1 )  ,TFMC  (3  ,  10  )  , TFM  (^  ,  100  ) 
COMMON/V/fiRCnv/VAPCO(c  ,1)  ,V*RC(7  ,10)  ,VAP  (3,  100  ) 
COMMON /W&LCDV/WALCnC 3  ,1 )  ,WAL ( 1 ,  100  ) 
7=0.0 

DO    3    I=1,10CC 
00     1    J=l,2 

ET( J.I  )=7 

1  CONTINUE 
n0    2    J=l,5 
EO! J,I)=2 

2  CONTINUF 

3  C.ONTINUF 

On  fl  1=1,100 
FOR(  1,1  )  =  7 
WA|  (  1  ,T  )  =  7 
00  4  J=l,2 
CTM  J,T  )  =  7 
PRM( J, I  )=7 
A     C°MTTNIJF 

0^    5    J=l,3 

FRFf J, I  )  =  7 

VAD( J,  T  )=7 
S    CONTINUE 

DO    6    J=1,A 

HIS( J,T  )=7 

TEM( J, I )=Z 
f    COMTINUF 

00     7    J=l,8 

NON( J,T )=7 
7    CONTINUE 
«    CONTINUF 

0"i    1*    1  =  1,10 

00    <?     j=it2 

EROC( J,! )=7 

NON'C(  J,  I  )=7 
Q    CONTINUE 

DO    10    J=l,3 

H!SC<  J,  I  )=7 

Tt:Mr(  j,t)=7 

10  CHMTJNUP 

DO    11    J=l,* 

cr>nr.(  j,i  )=7 

11  CONTINUF 

DO  12  )=1,7 
INTC  J, I )=7 
PRMC( J,I)=7 
V'PC( J,T )=7 

12  C^NTINUr 

DO  13  J=1,R 
ET*C( J,I)=7 

13  CONTINUE 

14  CONTINUE 

DO    le    J  =  l,3 

WALCO( J,l)=7 
1*    CONTINUF 

DO  20  J=l,6 

HI5C0( J,l)=7 
20  CONTINUE 

DO  16  J  =  1.5 

POGCO( J,l)=7 

FRRCO( J,l)=7 

INTCC( J,l)=7 

TEMCO( J, I )=Z 

VAPCP( J,l)=7 
16  CONTINUE 

DO  21  J=l,7 
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PRPCO( J,1)=Z 
21  CONTINUE 

TO  18  J=l,8 
BUMCO( J,l)=? 

pRMcnc  j,n=z 

13  CONTINL'F 

on  1<=  J=l,<? 

ETACO( J,l)-Z 
1<=  CONTINUE 

RETURN 

END 

SUBROUTINE    SFTUP 
REAL    MU 

DIMENSION  HI  ST 1(278)  ,HI ST2 (51 ) , HIST3 (375 ) , H IST4(  100  ), V 
1M/VT(  l5$) 

CnMM0N/EPGC0M/PEM(50)  ,PTE(50> ,DUES ,UE(50 , 50 ) , PE ( 50,  50 ) 
1,PHP.MU,T.DSIP(50),IPSIP 

C0MM0N/HISC0M/TZ,T1,T2,TF(25) ,C 1.C2 , C3. CA.BET A, ZM ( 5, 25 
1),XI( 50) ,HF(25,5) ,HUE ,HHUE , DLX2 ,C3M (50 ) , BET AM (50  ) 

C0MM0N/INTC0V/KP(2  0) ,KI N,KOUT ,MAT1 I , MAT  1  J, NET  A, I,IS,NS 
1, 1 T, N'T  I  MP ,NAM,NLEQ,NNLEO,NPNL,TTS,CASE(15),B(8),MWE,NO 
2N1ITE^,NITEMfKR17jNBT,NBT2tI0ENT,UMFStEREQ1EPS 

COMMON/PRMCOV/TIMF(5C) , PRE (50) ,  PTET  ,S  (  50  )  ,  ROK  AP  (  50  ) ,  RN 
10SE,VK4F,NDISC,IDISC(50) ,NSD(10) fMS Q( 10 ) , I  PR E.RADNO, CO 
2NEt0ELT, CYCLE 

COMMON /VARCCM/F (5,25)  ,FTM1 (5 ,25 )  ,FTM2 ( 5 , 25 ) 

C0MMON/WALCOM/FW(50)  ,RH0VW(50) , IFW, I RHOVW , JRHOVW 

EQUIVALENCE  (HI  SU  ,XI)  ,  (HI  ST2,PTF)  ,  (HI  *T3,P),  (HIST4,FW 
1),( VMAT,T7) 

2  FORMAT( 1H18X4HTIMFE12.5,56H    SECONDS    - 

I __.    -,6X,A5,3X,A2, 1H.,A2) 

3  F0RMAT(1H1»3X4HTIME,E12.5,34HSEC0NDS    *    *    *    STRE4MWISE 
1DIMENSICN    ,    E12.5,13H    FEET    ****/) 

A    FPPMfT(  1H18X4HCASEF12.5,56H -    - 

I ,6X,A5,3X,A2,1H.,A2) 

5    FnRHAT(lHlf2X12H*    *    *    *    CASE    ,E12.5,27H    *    *    *    STREAMWI 
1SE    DIMENSION    ,E12.*,13H    FEET    ****/) 
7771    FOPMAT( 1D8E1C.3) 
9001    F0RMAT(I3,7E10.3) 
J=MOD(  ITPM, 2)4-1 
IF(MWE)1C1,154,154 
101    CALL    REOASE 
NBT=3 
NBT2=o 

IS  =  1 
IT=1 
IF(KPCD)     104,  10*, 103 

103  CALL    LINMAT 

104  KR17=KR(17) 

15*    IF( IS+ITFM-2)    153,153,1572 
1572    IC(NITEM-1)1574, 1262, 1574 
1574    WRITF(NBT)HIST1.HIST2,HIST3,HIST4,VMAT 
IFfITCM-1)     157,156,157 

156  IDUM=NBT 
NBT=NBT2 
NBT2*IDUM 
REWIND    NBT 
REWIND    NBT2 
no    TO    155 

157  IF( IS-1)    153,153,155 

1*5    PEAD(NBT2)HIST1,HIST2,HIST3,HIST4,VMAT 

GO    TP     1262 
153    CONTTNUF 

IF(IT-l)     105.105,106 
C  INITIAL    GUESSES    FOP    PRINCIPAL    DEPENDENT    VARIABLES. 

C  CALCULATF     (KR(2)=0),     INPUT     (KR(2)=1),    OR    US F    VALUES 

C  PROM    PREVIOUS   CASE     (KR(2)=2>.    NOTE:     LATTER    REQUIRES 

C  SAME    ETA    VALUE.     ITS    UTILITY    IS    FOR    REPEATED    SIMILARITY 

C  SOLUTIONS.    IT    OBVIOUSLY    CANNOT    BE   USED    FOR    FIRST    CASE. 

105  CALL    FIR^TG 
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if6    CALL    RCFCPN 

KR( 2)=2 

IS=1 
1262    IP(TIMF(  I)  )1C53,105^,105^ 
10C3    TIMF( 1)=-TTME(1J 

WRTTF(K0UTf5)TTME( ITFV) ,S(TS) 

TIME< 1)=-TIME(1) 

GO    TO    126 
10c4    WP!TF(KOUTt3)TIW!P(  TTM)  f  S<  IS) 

W"TTE(6f3)TTyF(  !TEW)  ,S(IS) 

COMPUTE    HISTORIC     INFORMATION 
126    CALL    HI^TXI 

mwf=C 

START    OF    I^RA^ION    LOC° 
158    IT<:  =  0 

KR(  17)=KR1"' 
lc9    RCTUPN 

ENO 

SUBPPUTTVE    RFC;  CF 

COMMON /FPGC ON /PFM(^0)  ,pt«=(-"))  ,PUFC  ,UE(50,5Q  U°E<50te50) 
l.P.HO,  ^l)tT,nciP(c;o)  ,irST  D 

COMMON /FT ACCV/FTi(2c)  ,PFTA(2M TnSO(2M ,  DCUC  24), Rl( 2*  )» 
1B2<  2M  ,LAR(100)  ,BM(73,2e) 

COMMON/ I NTrrv/KP  (2r)  ,KTN,KCUTf  MjiTH  t MAT  1  J f NET  A,  I »  IS , NS 
1,TT,MTTMF  ,N» M,W|  POf NNLFO,NRNLf TTS»C 4-E( lc ) , B ( 3 ) , MWF t NO 
2Nf  ITFMfMTEMfKRl"',NeTfNBT2  1 1 DENT  ,UMFS  ,  FREO,  EPS 

COMMON  /PR  MCPM/t  IMF  (cO)  ,°RF(50)  ,PTFTtc  (50  )  ,c>0KAD(50),  RN 
10SFfVK*P,NDT?r.  f  IOISCCJO)  tNSo(10)  ,MS0(1C)  ,  I  PP  E  ,p  APNO,  CO 
2NE,DFL"r,CYCLF. 

COMMON /PR  PC CM /PP ,Sr ,XV(5) 

COMMON /W£LCOM/FW(cO)  ?°WPVW(50) ,  *FW  f  I  RHOVW  ,  J  PHOVW 

1  F3P.MJT<lHl///////l«?X<-°'-'ilNST!?ADY    BOUNDARY    LAYER    INTFGPA 
1L    MATQTx    PROCEDURE     // 

127X25Hnrc?RTMFNT  n«=  f  cooNAilTICS  /27X25HN^VAL  POSTGPAnil 
2&tc  c-mppl  /30X2CKMPMT^REY,  '"ftLTFOPMIA  // 32X  16^  I  NVFC  T  I 
3GATI0N    OF     //10X15A*, //3?X15 

4HC0NTRPL  NUMBERS  /11X5FH1  2  3  4  5  6  7  3  9  10  1 
51  12  13  14  lc  16  17  18  19  2C  /° X20 1 3// 35X 10HET A  VALUES 
6    /4( 5X1F7E10.3/) /) 

2  FORMATS 2CX17HNOSF     RADIUS  =    C10.5,13H    FT  ) 

3  FORMA T( 20X17HCONE  HALF  ANGLE  =  p10.5t13H  DFG  ) 
±  P0PMAT(20X17HDFMCITY  =  C10.S,13H  LB/FT  CU  ) 
c    cncMjT(  2CX17HVI  <:cnqTY                     =    1CE10.3,13H    LB    SEC/PT    S 

10     ) 
a    PHPMfTf 2CX17HTI wf     GRID  =    1"C10.3,13H    c EC 

1       ) 
7    FOP«AT(  2CX17MC  YCLF     TI  ^  =    1PF10.3,13H    «  FC 

1       ) 
o    FORMAT(  1W1//////////31  X1PHSTATICN    TNPIJT     PA"""  a     //1X11H0I 

I STANCE, FT       ,  12 X5 ( 1 XI ^E 1 0. 3 ) /^ (24X* ( 1X1 PE10 . 3 ) / ) ) 
9    FORMAT( 1X5HRCKAP     ,  IB Xc ( 1 XI P F10 . 3 ) /9 ( 24X5 ( IX 1 °E10 . 3 ) / )  ) 

10  C0PMtT( 1X20HPRFSCURF     PATIO,    MEAN    , 3X l*  ( 24X5 (  1X1 °F  10 . 3 ) 
15( 1X1PF1C.3)/) ) 

11  FORMAT(  /cx68HTEMPFPiTlJRE  TOTAL  MPAN 

1  ocrit.  GSCIL.     /21 X53HPPES SURF  VELO 

2CITY  AMFLITUDF    FREQUENCY    /7X66H(0FG    R)  (A 

3TM)  (PT/SFC)  PARAMETER  (HERTZ)    /5(5X 

A1PF10.3)///) 

12  FOPMAT(  1X21HPNTP0PY    PROP,BTU/LB    R     , 2 ) /Q ( 24X5 ( IX 1 PE  10 . 3 
1X5( 1X1PF10.3)/) ) 

13  PORMAT(  1X22HVASS    FlUXtLB/SEC    FT  **2     ,  1XC  (  IX  1PF  10  .  3  )  /<=  f  2 
14X5( IX1PF1Q.3) /) ) 

TF    (  ITEM-D3C2, 302,1^4 
3C2   CALL    INPUT 

WRITE (KCUT,1)CA SE , KP , ( F TA ( JN)  ,JN=1 ,NETA) 

WPITE(K0lJT,ll)Tf  PTFT,IJMFS,EPS,PREO 

NTTFM=NTIMF 

MTIMF=1 

COMPUTE  INFORMATION  MFcDFO  tq  CCNSIOER  DISCONTINUITIES 
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30  3    J-l 

102  "SP( i)=i 

inisci  n=o 

S(  1)=ABS(  san 
IF    (NS-1)     105.105,1021 
1021    DO     101    I<=2,NS 

misc(  is)=o 

1F(S(IS))     103,101,101 

103  NSD(  J)  =  IS-M<sD(  J)-H 
S( !S)=-S( IS) 
IDISC( IS)=1 

J=J*1 
^sn(j)=is 

101    CONTINUE 

105    NSD< J)=NS-MSO( J)  +  l 

ND!SC=J-1 

IP( TAB S(KR(6)-2) -1)207,208,207 
207    IF(NS-1)234, 234, 20^1 
20^1    IIS=2 

LN7  =  i 

TF(RCKAP(  1) )    223,226,226 

223  RAONO=-POKAP( 1) 
ROKAPC 1)=0. 

OH    229    IS=2,NS 

TF(ROK/SP(  IS)  )    22^,224,225 

224  IP(KR(6))    221,221,222 

222  ROKAP<  IS)=SUS)  *SIN(PADN0/57.2O573> 

on  to  22Q 
221  ROKAP< IS) =RAONn*SI N( S (I S) /RAONO) 
229  CONTINUE 

GH  TC  234 

225  IF  (IS-NS)  2251,234,234 
2251  IIS=IS+1 

LN7=IS 
22h    00233  IC=IIS,N^ 

IP(PPK&r( IS)  )  233,233,227 
227  IFUS-1-LN7)  232,232,228 
225  LN7=LN7+1 

ROKAP(LNZ)=ROKAP(LN7-l)  +  (S(LNZ)-<:(LN7-l)  )/  (  S  (  IS  )  -5  (LN7 
1-1) )*(ROKAP( IS)-RPKAP(LN7-1) ) 

GO  TO  227 

232  I N7=I S 

233  CONTINUE 

234  VKAP=1. 
GO    TC    210 

205    00    209    IS=1,NS 
2QQ    RGKAP( IS)=1. 

VKAP=0. 
21C    CONTINUE 
191    STEF    =    .481E-12 

IF(KP(6)-l)182,182tl83 
152    WRITE(KrUT,2)PtOMO 

WRTTP(KrUT,3)C0NE 

183  WRITE( KCUT,4)RH0 
WRITE(KCUT,5)MU 
WRITE(KrUT,A)DELT 
WRITF(KCUT,7)CYCLE 

WRITE (KOUT,P)f SCI) ,I=1,NS) 
WRITE(KCUT,9)(R0KAP(I ) ,1=1 ,NS) 
WRITE(KCUT,1C) (PPF (I) ,I=1,NS) 
WRTTE(KrUT,12)(D*IP(I) ,I=1,NS) 
WRITE(KCUT,13)(RHPVW(I) ,I=l,NS) 

184  RETUPN 
ENO 

SUBROUTINE    INPUT 

"PAL    MU 

COMMON /EOGCOV/PEM( 50) ,PTE(50) ,OUFS ,UF ( *0 ,50 ) , P* < 50, 50 ) 
l,RH0,MUtT.0SIP(50) ,IDSIP 

COMMON/ETACOM/ETA(25)  ,DETA ( 24 ) , DSO( 24) , DCU ( 24  ), 8  1( 24 ) , 
18  2(24)  ,LAR( 100)  ,BA 1(73, 28) 
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rPM*0N/INTC0M/KR(20l  ,  KT  N, KOIJT,  MATl  I  ,MATU?NFTA,  I,  IS,NS 
1,TT,NTIME,NAV,NLEQ,NNLFQ,NRNL,ITS,CA^E(15 ) , B ( 8) , MWF, NO 
2N,TTFMfNTTFMfKP17fNBTtNBT2  , I  OENT  ,IJMFS  .  PRFQ,  F°S 

COMMON /FRMCOM/TIMF  (50)  ,DRF(50) ,  PTET,C  (  50  >  .  R.OK  &°  (  50  )  ,  DN 
IPSE ,VK*F,NO! ~C ,1^1 SC(5C) ,NSn(10 ) ,M$D( 1C ) , I PRC , R A ONO, CO 
2NF,PFLT,CYCLfr 

CnMMON/VARCO/F<5,25)  ,FTMi  (5  ,25  )  ,ETM2  (  c  ,  25  ) 

COMMON  /WALCOM/FW( 5  0)  ,RH0VW(e0) ,  IFW, I RHOVW , JRHOVW 

1  FORMAT(20I1,15A4) 

2  FPPMAT(2I3,E10.M 

3  FORMAT*  8F10.4) 

4  FORM£T( 5(2XE14.7)) 
RP^O    (KIN, 1)KR, CASE 
P.FAD    (KIN,2)NTTMF 
READ    (KIN,2)NS 

RF£0    ( KIN, 3) (S( Ic) ,IC=1 ,NS) 

!p  (KQ( 1)  ) 103,103,100 
100  IF(KP(  12))101,1C1,102 
Fl    REAO    (KINT?)NETA 

Rr/0     (  KIN,  3)  (ETA  (I  )  ,1=1  ,NFT4) 

102  C&LL    STPETA 

1C3  IP    (KR( 6)-l)  10^,104,10* 

104  READ    (  K!N,3)rONF.  ,PMO^E 

1C5  IF( IABS<KR(6)-2)-lUG6,107,lG6 

1C6  REAO    (KIN, 3) (ROKAP(Ic)  ,IS=1,NS) 

IC  REAO    (KIN,3)PTFT 

1C71  I«=    (KP(2)-1)U2. 110,112 

110  REAO    (KIN, 4)  ((F(I  ,J)  ,J  =  1  ,NPTU  ,1=1,4) 

112  IFdTEV- 1)115, 113,115 

113  IPCTP=1 

TDDF=1 

00    114    1  =  1, NS 

114  DSIPf I )=0. 

115  IF(KR( 5)-5)  118,116  ,118 

116  ip(  rnsTP-iTCM)iiP,ii7fnq 
\\-f  oftr,   (  kin, 3)  mqp(! )  ,1=1  ,NS) 

pc»0     (  KTN,  ?)  Tr^I  P 

118  IF    ( IDPF-ITPM)12C,11Q,12  0 

119  READ    (KIN, 4)  (PPFU  )  ,T  =1,MS) 
RE^D    (KIN,2)IPPE 

120  JRHOVW=0 

IF(  ITCM-l)  121,121,122 

121  IPW=1 
tohOVW=1 

122  REAP    (KIN,3)T 

12c    IFf ITEV-1) 130,130,131 

130  IF(KP( Q)-l) 131,135,138 

131  IF( IFW-TTEM) 132,133,132 

132  IF( TFW-1) 134, 13*. 134 

133  REAO  (KIN, 3)  (FW(I  )  ,1=1 , NS ) 
REAO  (KTM,2)IFW 

134  GO  TO  13« 

135  IF(  IRHOVW-ITEM) 136  ,137,136 

136  GO  to  139 

137  RFA.0     (KIN, 3)     (RHOVWd  )  ,1=1  ,NS) 
REAO    ( KIN,2)IRH0VW 

138  PEA0(KIN,3)UMF^,FPE0,CRS 
JRHOVW=C 

IF    (FREO) 141,1^1,13° 
13c   CYCLE  =  1./FRFQ 

BAG=FLOAT(NTIME) 

OELT=CYCLE/*AG 

TIMF( 1)=0. 

no    140    I=2,NTIMF 
1A0    TIME(  I)=TTMP(I-l)*nELT 

GO    to    142 
141    timp( l)=C. 

OELT=Q. 

CYCLE=0. 
1^2    CONTINUE 

IF(KR(  16))  143,1^5,143 
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143    DO    144    1=2, NS 

1A4    P9P<  I)=PPE(I)*PHO*tJMFS*UMFS/136175. 57+1 
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RETUPN 
END 


SUPRHU 
COMMON 
1B2(24) 
COMMON 
NIT, NT 
2N,ITFM 
IF(KR( 
IF(KR( 
IF(KP( 
IM*R( 
NFTA  =  8 
E  TA  (  1 ) 
ETA(2) 
FTM3) 
ETA (4) 
ETA(5) 
FTA(6) 
E  TA  (  7  ) 
ETA( 8) 
GO  TO 
NETA=8 
E  TA i  1 
PTA( 
ETA( 
ETA( 
ETA( 
ETA( 
ETA( 
PTA( 

GO  TC 
NETA=1 
ETA(  1 
ETA(  2 
ETA(  3 
ETM  4 
FT/.  (  5 
ETA(  6 
ETA(  7 
ETA(  8 
ETA  (  «? 
FTA( 10 
GO  TG 
NFTA=1 
ETA(  1 
ETA( 
FTA( 
ETA( 
ETA( 
FTA( 
ETA< 
FTA( 
ETA( 
FTA(  10 
ETA  {  11 
ETA(  12 
GO  TO 
NFTA=1 
FTA(  1 
ETA( 
ETA( 
ETA( 
ETA( 
FTA( 
FTA( 
ETA( 
ETAC 


TINE    STDETA 

/FTACP«/ETM25)  ,DETA  (24  )  ,  DSQ(  24) ,  OCU  (  24  )f  R  I  (  24  ) , 

,LAR(1C0)  ,RA1(73,28) 

/INTCOM/KR<20)  ,KtN,KOUT,MATlITMATU?NETAf  I,  IS,N$ 

IMF  fNA  V,NLEQ,NNIEQ,NRNL,ITS,CASE(15),B(8)  ,MWE,NO 

,NITFMtKR17,NBT,NRT2,IOENT,UMFS,FRFO,EPS 

121-2)4,5,1 

12)-4)6,7f2 

12)-6) 8,9,3 

12)-8)10,11,12 

=  C.C 
=  C5 
=  1.C 
=  1.5 
=  2.0 
=  3.0 
=  5.0 
=  6.0 
13 

)=0.C0 

)=C30 

)=1.60 

)=2.40 

)=3.20 

)=4.00 

)=5.00 

)=6.0 

13 

0 

1*0.00 

)=0.60 

)=1.2C 

)=1.80 

)=2.40 

)=3.C0 

)=3.60 

)=4.20 

)=5.C0 

)=6.00 

13 

2 

)=0.00 

)=C50 

Ml. 00 

1=1.50 

)=2.00 

)=2.5C 

)=3.00 

)=3.50 


)=4 
)=4 


00 

50 


)=5.00 
)=6.00 
13 
4 

)'  =  0.00 


)=0.40 
)=0.80 


2 
3 
4)=1.20 

5 

7 


)=1.60 
=  2.00 
)=2.40 
)=2.80 
)=3.20 
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10 


1 1 


12 


FT/  ( 
r-  T/  ( 
FTA  ( 

ETA( 
ETA( 
FTA< 

F  TA  ( 


FTA(  10)=3.60 
rTA(  U)=4.00 

ETA ( 121=^.50 

FTM 13)=5.C0 
FTA(  14)=6.0C 

on  tc   13 

MCTJz  i«5 
FTM l)=C.OO 
ETA (2) =0.20 
FTM  3)=C.4C 
ETA<* )=C.60 

eta( 5)=r. ec 

ETA ( A)=1.CC 
PTA( 7)=1.20 
FTM  8)  =  1.40 
ETA(«)=1.60 
ETA( 10) =2. 00 
CTM  11)=2.50 
CTM 12)=3.25 
FT/.  (  13)=^.  CO 
Etm  14)=5.00 
CT£ ( 15)=6.C0 
GO    to     13 
N^T£=1Q 
ETA(     1)=C00 
CTM     2)=C. 3C 
3) =0.60 

5)=1.20 
6)=1.50 
7) =1.80 
8)=2.10 
S)=2.4T 
ETf ( 10)=2.7C 
FTM  H)=3.03 
FTM  12)  =  3.30 
FTM 13)=3.60 
ETA( 14)=3. CC 
ETA.(  lc  )=4.20 
ETA( 16)=4.60 
ETA(  17)  =  ?. CO 
-TM  18)  =  6.00 
GO  TO  13  • 
NFT£=22 
ETA<  1 )=C.000 
ETA(2)=C125 
FTM3)=C.25C 

F^Af  ?)=C, 

ETA (  6)=C 
ETA(7)=C, 
ETA( °)=C 
ETA( ~)  =  1 
ETA( 10)=1.125 
ETA(  11)  =  1.2C?C 
FTM  12)  =  1.3-?c 
FTA(  13)  =  1.  5CC 
FTM  14)=1.75C 
FTM  15)=2.00C 
FTA(  16)  =2. 250 
pr& {  17)  =  2.5CC 
FTM  18)  =  3.C0C 
ETAf 1 -)=3. 50C 
ETA( 20)=4.0CC 
ETA( 21)=C.0CC 
FTA{ 22)=*.Q0C 
GO    TO     13 
NETA=25 
ETA(     11=0.00 
FTM     2)=C.2C 
ctm     3)  =  C.4C 
FTA<     4)=C.60 


c00 
625 
7  50 
P7=; 

000 
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ETA(  5)=C80 
ETA(  6)=1.00 
ETA(  7)=1.20 
ETA(  8)=1.40 
FTA(  <3)  =  1.60 
FTA( 10)=1.80 
FTA(  ll)  =  2.00 
FTA( 12)=2.20 
FTA( 13)=2.4C 
^TA(  14)=2.60 
FTM  15)=2.80 
ETA(  16)=3.0C 
ETA( 17)=3.20 
ETA(18)=3.A0 
ETA( 19)=3.60 
ETA(2C)=3.8C 
ETA(21)=4.00 
FTA(22)=4.25 
ETA(23)=4.50 
FTA( 24)=5.C0 
CTA(25)=6.C0 
13  RETURN 
END 

SUBROUTINE    LTNMAT 

COMMON /ETACQM/ETA (25)  ,DETA (24 ) , DSQ( 24 ) , DCU { 24 ), Bl ( 24 ) ♦ 
1B2(24)  ,LAR( 100) ,RA 1(73,28) 

COMMON/INTCOM/KR(2  0)  , KI N  , KQUT ,MAT1 1 ,  MAT 1J. NET  A, I , IS ,NS 
l,IT,NTTME,NAK,NLEG,NNLEQ,NRNl,TTS,CASE(15),B(8),MWE,NO 
2N,ITEM,NITEMfKR17,N8T,NBT2,IDENTtUMFS,FREQ,E&S 

COMMON /TEMC0M/SP0UM2 (10)  ,DER(50 ) .DUMM1 < 25 ) .SLOPE (25  It R 
lE0UM(25)fSDUMl  (50) ,SDUM2(50) ,FWDUM(50) ,X ICON (50 ) , FWCO 
2N( 50),DUDS(5C) 

DO    104    I=2,NETA 

DETA( I-1)=ETA(I J-ETA(I-l) 

DSO( I-1)=DETA(I-1)*DETA(I-1) 

Bl(  I-1)=B(  3)*DS0(I-1) 

B2( I-1)=2.*B1(I-1) 
10 A    DCU( I-1)=DFTA(I-1)*DS0(I-1) 

MATlI  =  3*NETA-2 

MATU  =  NETA>2 

DO    105    1  =  1, MATH 

DO    105    J=1,MAT1J 
1C5    BAH  1,J)=0. 

BAK  1,1)  =  1. 

BA1(MATH,2)=-1. 

DO    106    I=2,NFTA 

IF  (KR(10))  1041,1041,1042 
1042  I*  (I-NETA)  1044,1041,1044 
1CM    DUM1=B(5) 

DIJM2=B(4) 

DUM3=B(3) 

DUM4=B(2) 

OUM5=B( 1) 

DUM6=B( 1) 

GO    TO    1043 
104*    DUM1=0. 

DUM2=B(3) 

DUM3=0. 

DUM4=B( 1) 

DUM5=0, 

DUM6=1. 
10*3    BA1    (  I-1,H-1)=DUM2*DCU(I-1) 

BA1    (  I-l,H-2)=DUMl*DCU(I-l) 

J=I+NETA-1 

BA1    (  J,H-1)=DUM4*DSQ(I-1) 

BA1    ( J,I+2)=DUM3*DS0(I-1) 

J=J+NETA-1 

BA1    ( J,I+1)=DUM6*PETA(I-1) 
106    BA1    (  J.H-2)=DUM5*DETAU-1) 
3060    F0RMAT(2X1P12E10.3) 
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o0f4  FPRM*-M?X37HLINF«P  MATRIX  FOR  MOMENTUM  EOU  AT  I  CNS  ,  I  3,  2H 
l  X!3,27H,  BEFOPF  AND  AFTER  SOLUTION! 
IF(KP(15l-l)  9062,9062,^061 
Q061  WRITP(KPUTf9064l  MftTHfM4TlJ 

00  ?063    1=1, MATH 

=  06  3    WRITF(  KrUT,9C6C)  (PA1  (I  ,  J)  ,  J  =  l  ,MAT1J) 
^062    00    107    LI=lfMATlJ 
107    CALL    MftTnNE(RAl(l,LT )) 

IF    (KR(  151-1)9072, °072,°066 
=  066    On    ^067    T  =  1,WA"T1! 

9067    WRITF(K0UT,QC6C)  (BA1  (I  ,J)  ,.J=i  ,MAT1J) 
9072    NLF0  =  MA^H 

NNLEQ=M/ T1J 

NAM=NNLF0-1 

NRNL=1 
9=02    CPNTINUF 

RPTUPN 

END 

SUBRPUTIME    MATONF(F) 

OTMENSICN    E<  1) 

CPMMPN/FTACOM/ETA(?c)  ,DETA(24)  , Or0(24)  ,00.1(24),  PI  (24), 
1B2( 24)  ,L*R(100)  ,BA 1(73,28) 

C0MMPN/INTCPM/KR(20) , KI N ,K CUT , MATl I , MAT1J,NETA. I, IS, NS 
l,IT,NTIVFtN&vtK«LFO,NNLcO,NRNL  , T TS ,C ASE ( 15 ) , B ( 8 ) , MWE, NO 
2N,!TeMtNITcv,KRl-r,NBTfNRT2  ,  I  OENT  ,IJMFS  ,  FPEQ,  FPS 

I=MATH 

DO    201    J=3,NETA 

1=1-1 
2C1    E( I)=F( I)+^( 1  +  1  ) 

I=NETA 

M=NETA+NETA-1 

K  =  C 

no    202    J=2,NET* 

1  =  1  +  1 
M  =  M+1 
K  =  K  +  1 

20  2    E(I)  =  E(I-l)-Etn+DETA{K)*E(w|) 

203  E(  1)=-F(  l)+nETMl)*F(K  +  l)  +  OSQ(l  >/2.*E(  j  +  i) 

1  =  1  +  1 
K  =  K  +  1 

no    204    J=3,NCTA 
1  =  1  +  1 
K  =  K+1 
v  =  M+l 

204  ^(M)  =  E(M-i)-F(M)+DFTMv)*E(K)+DS0(M)/2.*F(!) 
RPTnRN 

SUBROUTINE    RERAYfN.CDtNNN.IS  ,ND) 

niMENSTCK'  0(Nn,l)  ,Sn(2Q)  ,C(Nn,l)  ,L(2Q)  ,S(28),LL(2B),LL 
1L( 2P) 

CnMMCN/INTC0M/KR(20)  ,KTN,KCUT 

DTRECT    INVERSION    PPOOEDURE       —      C    IS    RFPLACEn    BY    C**-l 

Nl-N+1 

on    11    I=1,N 

LLL( I )=I 
113    L( ! )=I 
11    CONTINUE 

IX=-1 

IFfIS+2)     111,109,111 
1C6    F0PM£T(11H       L(I ) ,I=1,I3,5X     (3013)) 

107  FnPM£T(15H       (  (C (I , J)  ,J  =  1 ,I3,12H) ,<n<J)  , J=l, 13,    6H),T=1 
1,13, 15H)       BEFORE    PERAY) 

108  F0PMAT(2X1P11E10.3/(12X1D10«=10.3)) 

109  WRITC(K0UTtlC7)     N?NNN,N 
WPITF(KCUT,1C6)N,(L(U  ,1=1  ,N) 
IX=0 

no    110    1=1, N 
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11C    WRITF(*OUT,108) (C < I , J) , J=l ,N) , ( 0 ( I , J ) , J=l , NNN ) 
111    IS=-1 

TRIANGULATE    MATRIX 

DO    15    1  =  1, N 

00    160    M=1,N 

160  S(M)=ABS(C(I  ,M) ) 
IF(IS)     18,16,16 

18    1 5=0 

GO    TO    12 

REOUCE    ROW    I    BY    PRECEEDINC    ROWS 

16  DO    17    J=2,I 
K=L( J-l) 
DIV=-C(I,K) 
IF(DIV)    161,17,161 

161  C( I,K)=0. 
D0162    M=1.N 
DIVC=DIV*C< J-l.M) 
S(M)=&VAX1(S(M) .ABS(DXVCI) 

162  C( I,M)=C( I,M)+DIVC 
IF(NNN)     17,17,163 

163  00  164  M=lfNNN 

164  0( I,M)=0( I ,M)+DIV*0( J-1,M) 

17  CONTINUE 

SEEK  MAXIMUM  PIVOT 

12  OIV=0. 

DO  13  JJ=I,N 
M=L( JJ) 

IF(ABS  (C( I,w) )-0IV)13,13,121 
121  DIV=ABS  f C  <  T  ,M) ) 
K  =  M 
J  =  JJ 

13  CONTINUF 

SD( I)=DIV/S(K) 

t( J)=l(I) 

L( I )=K 

IF{  SD(  D-l.F-8)     130,130,14 

130  C( I,K)=0. 
IF($D<I))     131,131,12 
SINGULAR    MATRIX    RETURN 

131  IS=-1 

WRITE(KCUT,151)     (II  ,L(II)  ,SDUI),II  =  1,I1 

RETURN 
1A    DIV=C( I,K) 

C(I,K)=1.0 

K=LLL( J) 

LLL(  J)=LLLU) 

LLL(I)=K 

LL(K)=I 

NO&MALIZE    ROW 

IF(NNN)    1A3, 143,141 
141    DO    142    J=1,NNN 
1*2    D( I,J)=D( I,J)/DIV 
143    DO    15    J=1.N 
15    C  (  I  ,  J  )  =C  (  I  ,  J  )  /D I  V 

IF(IX)     152,150,152 

151  F0PMATC24H    °IVOT    ROW/COL/RES.  RATIO    5 ( 14 , 1H/ 13, 1H/E9 .2, 
11H,)) 

l^C    WRITE(KCUT,151)     (I ,L U ) ,SD < I ) , I =1, N) 
DIAGONALIZE    MATRIX 

152  NM=N-1 

DO    20    1=1, NM 
K«LU  +  1I 
DO    20    J=1,I 
OIV=-C( J,K) 
IF(DIV)19,20,19 
1?    C( J,K)=C. 

IF(NNN)     191,191,192 

192  DO    1<33    M=l,NNN 

193  0( J,M)=D( J,M)+DIV*D(I+1,M) 
191    DO    201    M=i,N 

201    C< J,M)=C( J,M)+DIV*C(I+1,M) 
20   CONTINUF 
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INTERCHANGE  COLUMNS 
DH  30  II=1,N 
T  =  I T 

21  J  =  L(I  ) 
L( I  )  =  T 

IF( J- I  )22, 30, 22 

22  TM  Ic)25,23,25 
2^  nn  24  M=l  ,N 

S(M)=C(M,I) 

24  C(M,T  )=C(M,J) 

I  =  J 

r,n    to  21 

25  IPf  T*-.l)26,?P,2fc 
2*  nn  27  m=1,N 

27  Cf MtT  l=f( m,  J) 

T=J 

GO  Tfl  21 
29.    00  2°  V  =  1,N 

29  C(M,I)=S(H) 
IS=0 

30  rnNTjMlJP 

IN~fpch'mgf  pdws 
nn  4  0  II=1,N 
I  =  T I 

31  j=u.m 

LL( T  )  =  T 

Ic(  J-T  )  32  ,40  ,32 

32  IPC  1^)35,33,35 

33  nn  it    v=l,N 

S  (  M  )  =C  ( I  ,  M ) 

34  C( !,M)=f ( JTM) 

TFfNNN)  343,343,341 
3^1  on  3*?  M=l,Nf^N 

SO(M)=pf  T  ,M) 
34  2  0( !,M)=P( J,M) 
3*3  IS=T 

!=J 

GH  tP  31 

35  !P( IS-J)36,38,36 

36  nn  37  m=i,n 

37  C(I,V)=C( J,M) 

IF(NMN)  3^3,373,371 

371  DO  37?  V=1,NKN 

372  0( T,^)=r( J,M) 

373  T=J 

on  to  31 

3°  On  39  M  =  l  ,N 

39  C(  !,*M  =  Sf  M) 
IP(NMN)  3°3,393,391 

3^1  00  3°2  M=lfN\N 
3C2  n( I,M)=?n(M) 
3«=3  IS=0 

40  CONTINUF 

IP(IX)  411,409.^11 
*07  «=0RM/iT(15H   (<C(I  ,J)  ,J  =  1,I3,12M),  (D(JI,  J=l,  13,  6H),I=1 
1,13, 15H)   4FTF&  PFPAY   ) 

409  WR!TF(KCUT,4C7)  N,NNN,N 
DO  MO  I  =  1,N 

410  WR!TF(KOnTf LC8) (Cfl ,J)  ,J=1 ,N)  , ( 0 ( I , J ) , J=l ,NNN  ) 

411  PPTURN 
END 

SUBPOUTTMP    FIRSTG 

COMMON  /FTACOV/FTM  25)  ,0PTA(24)  ,DSO(24)  ,0CU(24),  BK2M, 
1B2(  24)  ,LAR(  100)  , BAH  73,28) 

COMMON/INTCOM/KR(2  0)  , KI N,KOUT , MAT1 1 , MAT  1J.NET  A. I, IS, MS 
l,TT,MTTME,NAM,NLE0,NNLE0,NRNL,ITS,CASE(15),Bf 8),MWF,NO 
2N,!TFM,NItEM,KPl7,NRT,NBT2  ,  I  DENT  .IJMFS  .  FRFO,  EPS 

COMMON/PPMCOM/TIMffSO) ,°RE ( 50 ) , PTET ,S ( 50 ) , ROK \v ( 50 ) , RN 
10^p,VK/»P,NDTSC  ,101  Sf  (SO)  ,NSD(1C)  , MS  D(10)  ,  I  °RF  ,  R  AONO,  CO 
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2NF,0El_T, CYCLE 
COMMON/VARCOM/F(5,25) ,FTM1  (  5  ,25  )  ,FTM2 (5 ,25 ) 
IMK«(2)-1)     112,152,152 
112    00    129    I=1,NETA 
F( 1,1 )=0. 
OUMl=ETA(I )-4. 
IF(DUMi)     115,116,116 

115  DHM1=-C.1*DUM1 
DUM2=0.25*ETA(I ) 
F(2,I)=DUM2 
^(3,T )=OUMl 

GO    Tp    114 

116  F(2,I 1=1.0 
F(3,I)=C. 

114    Ff  4,1 )=0. 
128    CONTINUE 

IFf KR( 20) )152,152,200 

200  DO    208    1  =  1, NET* 
IF(ETA( I  )-5.)201,201,202 

201  F( 1,I)=ETA(I ) 
F(2,T)=ETA(I  1/5. 

F<  3,I)».l*(5.-ETAn)) 
GO    to    203 

202  F( 1,1 )=5. 
<M2,I  )  =  1. 
F(3,I  )=0. 

203  IF(ETMI  )-2.)204,204,205 
2C4    P(/*,I  )=-.l*ETA(I  ) 

GO    to    20  » 

205  IF(FTA(  I  )-5.  1206,206,207 

206  P(4,I )=-.067*(5.-ETA(I )) 
GO    TO    20  9 

207  P(A,I )=0. 

208  CONTINUE 
P( 1,1)=0. 
P( 2,1)=C. 
F(2,NETA)=1. 
F(3,NETM=C. 
P(4,NETA)=0. 

152    CONTINUE 
9O01    F0RMAT(4(2X7E10.3/)) 
qqqc    POPMAT( 2X19HDEBUG    FIRST    GUESSES) 
IMKR(15))    9902,9903,9902 

9902  TGNTINUF 
WRTTE(K0UT,9904) 

WRTTE     (K0UT,9901) ( (F(I ,J) , J=l ,NETA) ,T=1,4) 

9903  CONTINUF 
DETURN 
END 


SUBROUTINE    PSFCON 

RE^L    MU 

C0MMCN/EDGC0M/DEM150) ,fTE( 
1,RHC,MU,T,DSIP(50)  ,IDSIP 

COMMON/ETACO^/ETA(25) ,DETA 
19  2(24) ,LAR( ICO) ,BA 1(73, 28) 

COMMON  /HI  SC0M/TZ,T1,T2,T«=( 
lltXK  50), HP  (25,  5)  ,HUE  ,HHUE 

COMMON/INTCOM/KR(20) ,KIN,K 
1,IT,NTIMF,NAM,NLEQ,NNLE0,N 
2N,ITEM,NITEM,KR17,NBT,NBT2 

C0MM0N/PRMC0M/TIMEC50) ,PRE 
10SE,VK/>P.NDISC,IDISC(50)  ,N 
2NE,DELT,CYCLE 

C0MM0N/TEMC0W/SPDUM2(10)  ,D 
1EDU^( 25),SDUM1  (50),SDUM2( 
2N(50) ,DUDS(50) 

CnMM0N/WALC0M/FW(50) ,RH0VW 

1  F0RMAT(1H1//////////29X21H 
1TIMF=     1PE10.3.4H    SEC    /) 

2  FOPMAT( /5X68HSTREAMWISE 


50) , DUES, UE( 50, 50), PEC  50,50) 

(24),0SO(24),DCU(24),BK24), 

25)  ,C1,C2.C3,C*,9ETA,ZM(5,25 

,DLX2 .C3M (50), BET AM (50) 

OUT,  MATH  .MAT  1J.NET  A,  I,  IS,NS 

RNL,ITS,CASE(15),B(8),MWE,NO 

.IDENT,UMFS,FREQ,EPS 

(50) ,PTET,S(50),ROKAP(50),RN 

SD(10) ,MSD(10),IPRE,RADNO,CO 

ER(50),DUMM1(25)1SL0PE(25),R 
50)  ,FWDUM(50),XICON(5O),FWCO 

(50)  , IFW,IRHCVW,JRHOVW 
STATION   COMPUTED    DATA    /30X5H 


XI 


BETA 


82 


1  EDGE  STATIC  /5Xf  9HDI  MENS  I  ON        (LB/ 

2  VELOCITY  PRESSURE  /7X6 
35H(FFET)  SEC)**?  (FT/SEC) 
4         (ATM)  /) 

3  cnRVAT(5(5XlPE10.3)) 
ooOi  cnRM«T(^(5XlDE10.3) ) 
lx32    IF( TTEM-i)605t1538,605 
1*38  IOSIP=l 
I^PF=1 

no  i*4o  1=1, ns 

lc*0  D<!P( I  )  =  C. 

1*3*  I^  (KR(6)-3)1535, 1541 ,1535 
1535  IP  (PRF( l))1541t1542,1541 
lc42  On  1543  1=2, NS 

IP(PPF(I>)  1544,1543,1544 
15*4  L=I 

GO    TO     154  5 
15*3    CONTINUE 
1545   RMOSE«SfU/SCRT(  1. -PRE (D) 

HO    1*46    1=2, L 
15Afc    ope(I-1)=1.-S(T-1) /RNCSE*SU-1> /PNOSE 
15M    no    1*31    T=lfNS 

°TF( I )=PTET 
1531    t>cM(  I  )=PTET*PRF(I) 
8002    °I  =  3.141K<; 

no    8    1=1, NS 

nn    a    II=1,NITEM 

nDOS=(PTFT-PEM(I} ) /PHC 

n^DS=l3f 1?*.57*DPDS 

UE(  I,II)=*QRT(PPnS) 

UE(  TtIT)*UEI  I,II)*U.  +  FPS*COS((FPEQ*~(I)/lJF(  1,1  m  + 
1FL0AT( II)*2.*PI /FLOAT (NITEM) ) ) 

»E<  I,II)  =  PTFT-(PHC*UE(I  ,11  )**2./136  175.57) 
8    CONTINUE 
2Cr    IF(KP(15))    9?03,9«?0^,9P03 
QQ03    CHMTiNUF 

WPTTF(KOtJTtoooi)UF(Ic,ITFM)  t  PTC  US  )  ,T  ,RHO,  Ml) 
oqQi.    r*nNTTNUr 

c    forces  runs  to  zero  fcr  flat  platf  test. 

Ic    (KQ(7))8if83,«l 

81  OH    P2    JS=2,NS 

82  UF( JS,ITEV)=UF( JS-1,ITPM) 
°3    CONTINUE 

C  END    nF    EDGE    PROPERTY    LOOP,    START    OF    RETA    ANO    XT 

C  CALCULATION. 

6CK     XI( 1)=C. 

J=N0ISC  +  1 

DO    111    I 1=1, J 

K=NTP( T!) 

M  =  M$n(  T  T  ) 
LI.  =K  +  M-1 

IF    (  IT-1)     6052, 6C52,    402 
60c2    IF    (KP(6)-1)    40C,40l,^C2 
C  AXI^YMFTPIC    RLUNT 

4CC    On    403    I=M,LL 

SnUM2(  I)=SU  )*S(T) 
IF( S( T)  )    403,403,4031 
4031    01JM1=R0KAP(I  )/S(I) 

XICON(  I  )=UE<  I  ,ITFM)/S(I)*PHQ/4.*MU*DUM1*PIH1 
4C3    SPUMK  I  )  =  $0UM2U  )*^nuv2U) 
01=4. 

8ETAM(  1)=0.5 
GO    TO    AC6 
r  ©LAMAR    PLUNT 

4C1    no    404    I=v,LL 
S^UM2( I)=S(I ) 
TP(^(I))    404,404,4041 
4  041    XICONf  I)=UE(I,ITEM)/SU  )*RH0/2.*MU 
404    S0U«1(  !)=S(I  )*S(I) 
01  =  ?. 
PFT&M( l)=l. 
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GO    TO    406 
402    IF    (KR(H-2)    408,407,408 
C  AXISYMETRIC     SHARP 

407    DO    40°    I=M.LL 

SDUM2( I)=S(I )*S(I) 

DUM1=R0KAP(I  )/S(I) 

XirON(  I )=PHO*UF (I ,ITFM)*MU*DUM1/3.*DUM1 

409  SOUMK  I)  =  S(I  >*SDUM2(I) 

Kill)  =XIC0N(1)*S<1)*S(1)*S<1) 

0T=3. 

IF    (II-l)    4051,4-051,406 
C  PLANAR    SHARP 

4C?    00    405    T=M,LL 
S0UM2(  i)=sm 

XICON(  I)=PHO*UF(I  ,ITEM)*MU*ROKAPU)*ROKAP(I) 
405    SnUMK  I  )  =  S(T  ) 
01  =  1. 

IF    (II-l)    4052,4052, 4C6 
A0C2    XI( 1)  =XIC0N(1)*S(1) 

4051    MM=M 

FW( 1)=0. 
*C6    CONTINUE 

CALL  SLOPQ(K,S(  M)  , UE (*, ITEM), DUDS (M),DER(M)  ) 
C      FOPCFS  DUDS  TO  7EP0  FOR  FLAT  PLATE  TEST. 
TP  (KR(7) ) 94,86,54 
84  DO  85  JS=1,NS 

Q5  ouns(JS)=o.o 

86    CONTINUE 

IP    (KR(6)-2)    4066,4062,4062 
4066    Ic    (M-l)    4061,4061,4062 

4061  IF    (PNOSE)    4063,4063,4064 
r               OUFS    COMOUTgn    BY    SLOPO 

4063  DUES=0Unn  1) 
GO    tc    &065 

C  DUES    Fopy    EFFECTIVE    NCSE    RADIUS    US T NG    NEWTONIAN    PLOW 

4064  DUPS=SOPT{ 2. /RH0*PE(1 ,1 TEM) *32 . 17405*2 1 16 . 2 ) / RNOSE 
^065    XICON( 1)=RH0*MU/(2.*VKAP+2.)*DUES 

4062  CALL    SLOPQC K ,SDUV1 ( M) , XICON(M) , DER ( M) ,X I ( M) ) 
XF    (LL-MM)     111,4101,4101 

41C1    DO    MO    L  =  MM,tL 

410  8PTAM(L)=2./QI*XI (L) /UE ( L, ITEM) *S ( L ) /S 0UM1 ( L ) *DUDS (L ) / 
1XIC0N(L) 

111    CHNTINUE 
9803    FO°MAT(8E10.4) 
C  CALCULATION    CF    C3    MATRIX 

DO    13*    1=1, NS 
IF    (KR(6)-1)     137,137,158 
13^    IP    ( 1-1)     139,139,158 
13<=    C3MU  )=-SQRT(RFTAM(I  )  /(  DUES*RHO*MU  )  ) 

GO    TO     13  8 
15 R    C3M( I )=-S0RT(2.*XI (I) ) / ( RHO*ROK AP( I )  *UE  ( I ,  IT  EM)*MU  ) 
139    CONTINUE 
joh0VW=1 
^043    IF(JPH0VW)7047, 7047, 7Q45 
C  CALCULATE    FW    IF    PHOVW    GIVEN 

70^5    J=NDISC+1 

DO    7046    11=1, J 
K=NSD( II ) 
m=msh( ii ) 
LL=K+M-1 
no    20C    I=M,LL 
IF    (KR(8))     7049,7049,2291 
7049    RHOVWf I )=RH0VW(I)*C3M(I) 
2291     IF    (II-l)     7048,7048,    230 
7048    IP    (KR(6))     229,229.230 
C  VALID    AT    AXISYMETPfC    STAGNATION    POINT    ONLY 

22°    FWCON( I)=-RHCVW(I) /C2.*C3M(I)  ) 
!P(     1-1)    20^,209,232 
C  MODIFICATION    FOR    AXISYMETRIC    BLUNT    *WAY    FROM 

C  STAGNATION    PCI  NT. 

232    FWCON(     I)=*=WCON(     I)/S(     I)*ROKAP(     I) 
GO    TO    209 
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VALID  FCP  ALL  PLANAR 
230  FWCOMC I )=-RHCVWU) /C3MU) 

JF(KR<6)-2)  209,236. 209 

MHPIFICATION    FOR    AXISYMETRIC    SHARP 
236    ^WCriNH     T)=nWC0N(    i)/S(     I)*ROKAP(    I)/2, 
2C<=    CONTINUE 

PWP!JV(  1)=FWCCN(1)*S<1  ) 

IF(KR(f)-2)    241,237,241 

MODIFICATION    FOP    AXISYMFTRIC    SHARP 
23^    ^WPU^f  1)=FWDUM(1)*S(1) 
2*1    COMTINUF 
7046    CALL    SLCP0(K,SDUM2(") ,FWCON(M) fOER(M),FWDUM(M)) 

OH    126    T    =1,NS 

IF( T    -1)    124,12^,123 
12*    IF    (KR(f)-l)     133,133,123 
133    IP    (  5f  I  ))    113,113,123 
113    FW(  1)=PHCVW( 1) 

GO    T0    126 
12  3    FW( I )=PWPUV( I)/S0PT(2.*XI( I) ) 
126    CONTINUE 
70^7    WRITFf KrUT,l)TIMF( TTFy) 

WRITE f KPUT.2) 

OH    ^050    1=1, NS 
70? C    WRITP( KnjTt3)S( I ) , XI ( I ) ,BET* M(l ),UE( I ,  I*<"EM) ,  PEU  ,  ITEM  ) 

OP"HJPN 

ENP 

SUBROUTINE    SLOPQ(N,X,Y,S,Z) 
PIMFN'SIPN    X(  1)  ,Y(1)  ,S(1)  ,7(1) 
TF(N-l)     c,g,S 
9    S(2)  =  f Y(2)-Y<1) ) /< X(2)-X(l)  ) 
~t  1)=S(2) 
OC=S( 2) 
DO    7    1=1, N 
TP( 1+1-^)2,1,6 

1  OR=OC 

IF     (T-2)7,6,5 

2  XOT=x( T)-X(I+1) 
XTT=X(  H-l)-X(T  +  2) 

xTr  =  x(  !  +  2)-xu ) 

A?=Y( I )/( XCT*XTC) 

XnTT=XPT*XT"r 

AB=Y( 1+1) /XOTT 

/*r  =  Y(  T  +  2)  /(XTT*XTP) 

AAA=ftA*XTT 

*RP=A3*XTfJ 

ACC=AC*XCT 

OA=OC 

QB  =  S( T ) 

oc  =  sr i+i) 

S(  I  )=A«*(  XTQ-XPT)*ABB-ACC 
S(  H-1)=AB*<  XCT-XTT)  +  *CC-AA4 
S( T+2)=*C*( XTT-XTO)+A4i-A8R 

3  IF(  1-2)7,5,4 

*    S(  !  )  =  <  <(I  )+OA)/2. 

5  S(I)  =  <$U  )+Q8)/2. 

6  XP=X(  I)-X(  1-1) 
YS=Y(  I  )+Y( !-l) 
SO=S(  !  )-S( 1-1) 
SS  =  M  I  ) 

Z(  I  )  =  7.  <  I-l)  +  XD/2.*(  YS-XD/6.*SD) 
S(I)=SS 

7  CnNTINUF 
9    RETURN 

ENP 

SUBRPUTINF  HISTXI 

REAL  ML) 

DIMENSION  FD(4) 

COMMON /FDGCOM/PFM( 50)  ,PTE(50) ,DUES ,UE ( 50 ,50 ) , PE ( 5  0, 50 ) 
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l.RHO,  MU,T,DSIP(50)  ,  I  OS  IP 

COMMON /EtACOM/FTA< 25) ,OETA  (  24  )  ,  OSQ(  24)  ,  DCU  (  24  ),  Bl  (  24) , 
1B2C24)  ,LAR(100),BA1(73,28) 

COMMON/HI SC0M/T2.T1,T2,TF( 25) , C 1 ,C2, C3.C4.BET A, ZM (5, 25 
1).XI<  50)  ,HF(25,5)  ,HUF ,HHUE ,OLX2 ,C3M<50 ) , BET  AM (*0 ) 

COMMON/INTCOM/KR(2  0)  ,  KI  N,K0UT  ,MAT1 1  ,MAT1J,NETA,  I,  IS,NS 
l,IT,NTTMF,NAM,NLEO,NNLEO,NRNL,ITS,CASE(15),B(8),MWE,NO 
2N,  I  TFM, MI  TcytKRl 7    NBT,NBT2f  I  OENT, IJMFS.FRF0.EPS 

COMMON/PRMCOM/TIME(50) .PRE(50),OTETfS(50).ROKAP(50),RN 
10SF,VKAPtNOISC,IDISC(50)  ,NSD(10)  »MS 0(101*  I  PRE, RAONO, CO 
2NF,OELTfCYCLE 

COMMON /VARC0V/F( 5,25)  ,FTM1 (5 ,25 ) ,FTM2 < 5 , 25 ) 

INITIALI7F    AXIAL    VARIATION    TERMS 

NUL  =  0 

OLX1=0. 

DO    120    1=1,4 

120  PO( I )  =  C. 

IP  (KR(6)-1)  400,400,401 

400  !c  ( IS-3)  154,153,153 

401  IF( IS-2)154,153,153 

lc3  IF(KRO)-l)  154,143,143 

143  DLX1=ALCG« XI(IS)/XI(IS-1)) 
Ic  (KR(6)-1)  402,402,403 

402  IP  (IS-3)  154, 155, 157 

403  !P( IS- 2)  154,155,157 
157  IPfKP(3)-l)  154,155,144 

154  07=0. 
01  =  0. 
02=0. 
M=NPTA-1 

00  140  1=1,4 

DO  140  J=1,M 

7M( I,J)=0.0 
140  CONTINUE 

DO  141  I=1,NETA 

DO  1*2  J=l,3 

HM  I,J)=C. 
142  CONTINUE 

HFf I,4)=C. 
1M  H«=(  I,5)=C. 

OLX2=0. 

GO  to  130 

COMPUTE    TWC-    OP    THREE-POINT    DIFFERENCE    RELATIONS 

155  07=2. /OLX1 
01=-D7 
02=0. 

GO    TO    145 

144  J=TDISC(IS-1) 
IF( J)121, 121.155 

121  DZ=DLX1+0LX2 
Dl=-DZ/(nLXl*DLX2)*2. 
D2=PLX1/(D7*0LX2)*2. 
DZ=-Dl-r2 

1^5    0LX2=0LX1 

122  F0( 3)=D1*F(4,1)+D2*HF(1,4) 
DO    147    I=2,NETA 

FO< 1)=D1*F(2,I )+02*HF(I ,2) 

FD(2)=01*e<3,I)+D2*HF(I ,3) 

P0(^)=F0(3) 

F0(3)=ni*F(4,I )+0  2*HF(I ,4) 
147    CALL    TAYLORfDETA(I-l)  ,FD(2  )  ,  FD<  1 ) ,  ZMU,  1-1 )  ) 

SAVE    HISTORIC    VALUES 
162    HO    164    I=1,NETA 

HP(I,4)=F(4,I) 

HFf I,5)=D1*F(1 ,1 )+D2*HF(I, I) 

00  164  J=l,3 

HF( I,J)=F( J, I) 
164  CONTINUE 

COMPUTF  GROUPINGS  WHICH  DEPEND  ON  OZ 
13C  Cl=i.-»-D7 

BETA=BETAM(IS) 

C2=-2.*(C1+BETA+DZ) 
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C3=C3V(  IF) 
C4=RHP*MU*C3*C3 

CnMPHTF    TWO    OR    THRFF     PPTNT    0IFFFReNCF    TIME    DFRIVATIVFS 

Ir(fR( M-1)20G,2>~>0,201 
2C0    IP( ITFM-2)202,203,2C4 
2C1    T^r Ttpm-2)202, 203,203 

202  T7=o. 
T1  =  0. 
T2=0. 

GO   to    2CC 

203  T7=1./PFLT 
Tl=-T7 
T2=0. 

nn    JO    ?CC 

204  t7=1.5/D£LT 
Tl=-2.  /^c!_t 
T2=.5/nFLT 

SAVE    H!cTpPTC    VALUFS. 
?CC    nP    2CC    I  =  1,NET/S 

TM I)=T1*FTM1( 1 f I)+T2*FTM2 (1 ,1) 
206    CONTTNl.T 
2C"7    CONTINUF 
?or*.    F^p,  ma  T('XT2,6XI2/(  8X1D10F10.3)  ) 

!F(Kp(17))     <?<;0Kf<?906tQO0e; 
02CC    CONTTNU? 

WRITE (*CUTt9507) 
gaQ7    cnPMfiT(2X2*HPEBnG    IFtPLXl     ...     7«,HF?TF) 
f.]F"Tft  Ml=^f  TA-1 

W»!TP(  Kr()T,qco/.)     I  S»!TFMtPLXl  ,0LX2,07f01f  D?,TZ,T1,T2» 
10.  ltC2,C3,04tFpt  (  (  zm(I  f  J)  ,  j=l  tNFT*Ml)  t 1=1,4),  (  (HFf  lt  j  )t 
2J=1,C)  ,I=1,NETA) 
=-06    CONTINUE 
RETURN 
ENO 

SUBPnUTTN'r    T£YL0P     (DfFV,F,9) 

DIMENSION    FM(l),F(l),o(i) 

rnMMON/INTCON/KO(po)  t  KI  N  fKCUT  ,  mati  j  ,MAT1J,NFTA,  I,  IS,N5 
1 , 1  T , MT  T  vc  rvj-s  m  ml  r o  ,  NNLF  0 , N^ NL  f  I  T 3  , C  A?  E  ( 1*  ) ,  B  (  8 )  ,  MW  F,  MQ 
2N,!TFW,MTnM,KP17,Kq^,N^T2f  IPENT,fUMcS,FCFO,FPS 

r»2=n*n 

!F    (ko(  1C)  )     1,1,2 

2  IF    (  I-MFTA)    4,1,4 

A   fp=o. 

PC  1)  =  (  f  (FM(3|  /*.-FD/24.  )  *P-F(2)  '2.  )*D+F(1)  )*0 

P(2)=( ( (fp  /3C.-FMJ  3)  /R.)*P+f<2) /3. J*n-F(l )'2.)*02 

PC  3)  =  C. 

o(4)=(((CM/2C.-Fn/72. )*D-F(2) /*.)*P+F( 1 )/6.)*D2*0-P<3) 

1(3) 
GP    tp    3 
1    cn=P( 3)_cm( 3) 

PC  1)  =  (  (  (M3)  /6.-FP/2*. )*n-F(2)/2.)*P+F(l) )*P 

PC  2)  =  (  (  («rP/3C.-F{3)/«.  )*p+f(2)/3.)*P-c<1)/2.)*02 

p(  *)=-((  CFP/252.-F  (3)  /72.)  *P+F(2)/30.)*O-F(  1)/2<V.)*P2* 

IP 
-(3)  =  (  (  (F(  3) /20.-FO/72. )*P-FC2) /9. )*0+f( 1  }76. )*D2*D-P< 

14) 

3  CPWTNUF 
RETURN 
EN|p 

SURRPUTINF    ITERAT 

common /pump,  ov/eijMp,coRVAf  5  as  f,  icorm,woot,  1777, isp,  ix 

CnMMPN/r:RPcnw/ELc(73)  , FLA  (73)  ,  F LFMt  ELM  (  3  )  ,  EL*",  IFLM, NE 
1LM,  TLMM,pFL(73)  ,FNLF  C2B)  ,ENL(23  )  t«=NL5M,  ^NLMC  3  ) ,  FNL^M,  I 
2FNLMt\FNLM,iNLMM.PFNL(28) ,PRNL(M 

0P^MON/!NT0nv/KP(2O)  , KI N ,K0UT .MA^l I . MAT 1J . NET  A. I, I«,NS 
1,IT,NTTM*,NA  v^LEO.NMLEO.N^NJLi^S^A^E*  lc  )  ,  B(  8)  ,  MWE,  NO 
2N,TTFyfN!TcM,KRi7  tNBT,NBT2  ,  I  DENT  ,UMF<;  f  FRFQ,  E°? 

r.nMMpN|/VARCnv/F(e;  ,25)  ,  FT*n  («5  ,2^)  ,FTM2  (  c  ,  25  ) 
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1    FORMAT(22H    NCN-CCNVFRGFNT    OUTPUT    /) 
5    FnRMAT(6Xl3,lXF9.3,2(lXF9.4).2(2Xl3f2XlPF10.3)) 
7    FORMAT* /10X62HITERATED    VALUES  DAMP  MAX    LI 

1NEAR  MAX    MOMENTUM    /6X62HITS  TIME  FPPW 

2  ERROR  ERROR) 

CLOCK8=0. 
181     IT$=ITS+1 

CALL    LINCER 
NON  =  0 

CALL    NONCER 
PPPW=F(3,1) 

IF(KR(4)+KR(  17)+KP(19)+N0N)189,189,192 
189    IF(ITS-1)192,192,194 
192    WRITE(KOUTf7) 
194    CLOCKl=C. 

RTIMP=CLOCKl-CLOCKB 
8T!ME=.01*RTTME 
CL0CK8=CL0CK1 

WPITE(KOUTf5)ITSfRTIME.FPPW,EASE,IFLM,ELMMf IFNLM,FNLEM 
1°21    IP(ELMM+FNLMM+.001/(EA§E*RUM°)-.00  2)162,162,159 
162    NON  =  0 

RETURN 
159    IC(ITS-5C)     161,160,160 
16C    WRTTE(K0UT,1) 
WRTTE(6,1) 

IF    (ELMM+FNLMM-.C2)     162,162,1601 
16C1    N0N=1 
RETURN 

ITERATE    OP    OUTPUT 
161    !F(KR(M)     181,181,193 
1C3    M0N=-1 
RETURN 
END 


SUBROUTINE    LINCER 

C0MM0N/ERRC0M/FLFC73)  ,ELA(73)  ,FLEM, ELM ( 3 ) , ELMM, I  ELM, NE 
1LM,  ILMM,DFL(73)  ,FMLE(28)  ,ENL(28)  ,FNLEM,  ENLM(  3  ) ,  ENLMM,  I 
2PNLM.NENLM.TNLMM,DFNL(28I ,DRNL(6) 

COMMON  /ETACO*/FTA ( 2C )  , DETA ( 24 ) , DSQ ( 24 ) , OCU ( 24 ) , B 1 ( 24 ) , 
18  2(24)  ,  LAP  (ICO)  ,  BAH  73, 28) 

C0MM0N/INTCC*/KR(20)  , KI N ,KOUT ,M in \ ,  MAT U . NET  A, I , IS,NS 
l,IT,NTIMF,NAM,MLEO,NNLEO,NRNL,ITS,CASE{15),B(8),MWEtNO 
2M,TTEM,NITEM,KR17fNBT,N8T2 . IDENT,UMFS , FREO, EPS 

COMMON /PR MCOW/TI ME (50)  ,PRE  (50  ) ,  PTET  ,S  (  50  )  ,  ROK  AP  (  50  ) ,  RN 
10SF,VKAPfN0ISC , IPISC(50) ,NSD( 10 ) ,MS D( 10 ) , IPRE,RADNO, CO 
2NF,DELT,CYCLE 

C0MM0N/TEMC0M/SPDUM2(10)  ,DER( 50 ) ,0UMM1 ( 25 ) .SLOPE ( 25  ) , R 
1EDUM( 25) ,S0UM1  (50) ,SDUM2(50) ,FWDUM(50 ) ,X ICON (50 ) , PWCO 
2N( 50) ,DUD$(CC) 

COMMON /VARCOW/F( 5,25) ,FTM1 (5 ,25 ) ,FTM2 ( 5 , 25 ) 

EVALUATE    LINEAR    ERRORS    FOR    MOMENTUM    AND    ENERGY 

NELM=2 

NENLM=NELM 

DO  401  T=2,NETA 

IF(KP(1CI)    4000,^000,4001 
4001     Ic    (I-NETA)    4003,4000,4003 
^000    DUM1=B(^) 

DUM2=8(5) 

DUM3=B(2) 

DUM4=B(3) 

DUM5=B( 1) 

DUM6=1,0 

GO    to    4C02 
4003    DUM1=B(3) 

DUM2=0. 

DUM3=8( 1) 

0U"4=0. 

DUM5=1.C 

DUM6=C. 
400  2    »=LP(  I-1)=-(F(1,I-1)+DETA(I-1)*F(2,  1-1 )  +  DS  Q(  1-1  )*B  ( 1  )*  F 
1(3,I-1)+DCU(I-1)*(DUM1*F(4,I-1)+DUM2*F(4,I) )-F(l,I) ) 
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M=I+MFT.A-1 

PLF(  V)=-<F (2.1-1 )*DFTA(I-1 )*F(3,!-l)*DS0< 1-1 ) *( DUM3*F ( 
14,!-1)+DUM4*£( J,I) )-F(?,I)  ) 
M=T-»-2*NETA-2 
401    FLMM)=-(F(3,I-1H-DETA(I-1)*0UM5*(F(*,  1-1  )+DHM6*  F(  4,  I) 
1  )-cf  3,1)) 
FLF1NETA)=-FC2,1) 
C  OF^PPMINF     MAXIMUM    LINEAR    ERRORS 

^C4    CALL    ABWAXfMATIT  ,FLF,ELMMt  !C|_H) 
C  FORM    oornUCT    OF    A**-l    ANr    LINFAR    ERRORS 

46<?    CU.L     MATPNF(FLF) 
473    CONTINUF 
RETURN 
END 

SUBROUTINE    ABMAX(N,X,XMf I) 
DIMENSION    X(l) 
1  =  1 

XM=ABS    (X(l)) 
Ic    (N-l)    4,4,5 
5    on    3    J  =  ?,N 
XT=AB^     (XM)) 
TP(XM-XT)     2,3,3 

2  XW=XT 
T  =  J 

3  CONTINUF 

4  XM=Xf ! ) 
RETURN 
E*'D 

SUBROUTINE    NCNCFR 

PEAL  MU 

OIMENCTTN  CORAP(25) 

COMMON /RUMCOM /R1JMP  ,C  OR  V4,E4SC,  I  CORM,WDOT,  P77f  IS*3,  IX 

Common /FOOCO"/PEM(^0) ,PTE<50) , DUES ,UF ( 50 ,50 ) ,©E( 50,^0) 
I^HH-mu    Tf0SIP(50)  .IDSIP 

COMMON  /ERRCO*/FLM 73)  ,FLM73)  ,P LFM, ELM ( 3 ) , EL MM, IFt_M,NE 
lLM,TtMMtnEL( 73) ,FNLFf28) ,ENL(2B ) ,FNL6M,ENLM(3), ENLMM, I 
2CNLM,NFNLM,INLMV,DFNL(2  8) ,0RNL(6) 

CnMMC\i/ETfiC0V/ETM?5>  ,OETA  (  24  )  ,  DSO(  24  )  ,OCU  (24),  Bl(  24), 
1B2(  24)  ,LAR(100) ,BA  1(73,2 8) 

CnMMON/H7SCOM/T7tTlfT2  ,tf(25)  ,C  1  ,C2  ,  C3  ,  C<<  ,  BET  A,  7M(5, 25 
1)  ,XT(  5C)  , HP  (25,5  )  ,HUF  ,HHUE,DLX2  ,C3M(50)  ,8ETAM(*0 ) 

COMMON /JK-TCCM/KR  (2  0)  ,  KT  N  ,KO'JT  ,  m^TI  I  f  MAT  1J,NET  *,  I,  IS,NS 
l,!T,MTIvr  ,NAVtML«r0,NMl.frQ,NRNL,TTS,CASE(15  ) ,  R  (  9  )  ,  MWE,  NO 
2N,ITCMtNTTCMfKR17,NRTf  N9T2  ,  I  DFNT  ,IJ"<FS  ,  FRFQ,  E°S 

COMMON/NONCOM//! "(2  8,2  8)  ,DVNL(2P  )  ,C5  ,  C6  ,C7  ,C8  ,  C9,  C  10,  C  1 
11,C12 

COMMQN/roMCOM/^IME (50) ,P0E(50) ,  PTET  ,S  (50  )  ,  ROK  AP  (  53  ) ,  RN 
10«F,VKAP,N0I«C , iriSC(5C) ,N*D(1C ) ,MS0(10) , I °RC , R AONO, CO 
2*'E,PELT.  CYCLE 

COMMON /PR PC Q* /PR, SC,XP (5) 

C0MMPM/TFMCGM/conuM2  (10),DER(50),0U*mi.(?e),SL0PF(25),R 
1cn,|M(  2?)  ,SOUMl  (EO)  ,S0UM2(50)  ,FWDUM(50)  ,XIC0N(5C  ),cWCO 
2N( c0) ,0UPS(5C) 

roMM0N/VtPCOV/F(5,25)  ,PTM1(S,25),CTM?(C,25) 

CnMMON/UALC:v/PW(50) ,RH0VW(50) , !FW, IRHOVW , JRHOVW 

EOUIVALFNCE     (C ORAR (1 ) , AM( 1 ) ) 

EA^F=AMTNHFASF*2.  ,1.  ) 

IF    ( T  Tc    -    i)     10,c,10 


AND  ERRORS  FOR  NONLINEAR 


K 

EASF=1. 

RUMP    =    l.C 

4 

IC0PM=1 

COPM4      =      l.C     4-     10 

10 

CONTINUE 
IX  =  0 

c 

FVALUATF    COEFFICIENTS 

C 

EQUATIONS. 

c 

INTTI*LI7E     4*    MATRIX 
Dn    IE     I=i,NNLEO 
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FNLE( I)=C. 

DO    iR    J=1,NNLEQ 
15    AM( I t J)    =    C. 
C  EWL.    GROUPING*    WHICH    CHANGE    DURING    ITERATION    BUT    ARE 

C  MOT    F(fTA ). 

IF( ITEM-2)30,3l,32 

30  r.^=o. 

GO    TO    33 

31  C5=TZ*ALCG(UE(I  S,ITEM)  )  +  T1  *ALOG  (UE  (  ISf  ITEM-H  ) 
GO    TO    33 

32  C5=TZ*AL0G(UE(  IS,ITFM)  )  +T1  *AL0G1 )  )  +T2*AL0G(U  E  ( IS  ,  ITEM- 
l(UE(IS,ITEW-2)) 

33  CONTINUE 

C EVAL    COEFS,    AM,    AND    ERRORS    FOR    EDGE    BOUN    CONDS 

C  NEXT,       EVALUATE    NONLINEAR    ERRORS 

45    FN! E(NETA*l)=-( l.-F(2,NETA) ) 
FNLE(NETA+2)*-F(3,NETA) 
FNLEC  1)=FW(IS)-F(1 ,1) 
C  NEXT,    EVALUATE    ORIGINAL    AM    COEFFICIENTS 

AM( 1,11=1. 
AM(NETA+2,2)=-l. 

FMLF(NFT/i  +  i)=FNLE(NFTA+l  )  +  FLE  (2*NET  A-i  ) 
DO    50    J=1,MAT1J 
*0    *M(NFTA+1,  J)=AM(NETA  +  1  ,J)*BA1  (2*NETA-1,J) 
C  START    OF    MAJOR    DC    LOOP    FOR    CVAL    nF    COEFFS     AND    ERRORS 

C  AT    PACH    STATION. 

DO    120    I=1,NETA 
C  TPST    TO    BY°ASS    COMMANDS    THAT    CANNOT    BE    PERFORMED    AT 

C  FTA(l). 

IF    (1-1)75,75,55 
c5    CALL     IMONE 

C FVAL    GROUPINGS    WHICH    ARE    USEO    AT    1-1    AS    WELL    AS    AT    I 

75    C6=HF( 1,5) 

C7=F(3,I)+F(2,I)*(Cl*F(lfI)+C6) 
CP=C4*(TZ-C5) 
Co=Cl*F(2,I>-C8 
C10=C6+C1*F(  1,1  ) 
C11=-CA*tf(I  ) 
C12=C9*F(  1,1  ) 
IF    (1-1)     100,105,100 
C  BACK    TO   CONSERVATION    EQUATIONS 

100    CALL    IONLY 
105    IF    (KP(17))     120,120.115 

110    F0PMAT(21H    ALL    THE    COEFFIC IENTS / (1X1P12E10. 3 > ) 
115    WR I TE ( KOUT, 1 10 )C1,C2,C3,C4, 05,06,07,08,0^,010,011,0 12 

IX=-2 
120    CONTINUE 

1*5    IF(KP(1<5))     170,190,170 
170    CONTINUF 

WRITE (KOUT, 1751 
175    FORMAT(2X10HDFBUG    FNLE    ) 
IPO    FORMAT(/( 12X1P1CE10.3)) 

WRITF(KOUT,lflO) (FNLE (J) , J=1,MAT1T) 
1Q0    CONTINUE 

DO     1»<?    111=1, MAT1J 
1P9    ENL( I  I T  )=FNLE(III) 

CALL    RERAY(MAT1J,AM,FNL,1,IX,2P ) 
IF    (KP(17))    245,265,245 
2^5    CONTINUF 

WRITEC<nUT,255) 
255    F0RMAT(2X9HDEBUG    FLE     ) 

WRTTE(K0UT,180) ( FLE ( J ) , J=l ,MAT1 I ) 
265    CONTINUE 
C  DETERMINE    MAXIMUM    NONLINEAR    ERRORS 

CALL    ABMAX(MAT1J,FNLE,FNLEM,IFNLM) 
ELMM    =    ABS(ELMM) 
ENLMM=ABS(FNLEM) 
C  EVALUATE    NONLINEAR    CORRECTIONS 

DO    615    I=1,MAT1J 
615    DVNK  I  )=ENL(  I) 
625    CONTINUE 
C  EVALUATE    LINEAR    CORRECTIONS 
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on  *30  1  =  1, math 

00    <>30    J-1,MAT1J 
630    «=t.  F  (  I  )    =    FLEd)     -    OVNL(J)     *    8A1U,J) 

F|_E(NFT* )=0. 

00    62C    T=1,NETA 

L*!+NET£-1 
62C    COP*R(  !)=FLE(L) 
6*5    CONTTNUF 

ir    {.33    *■    CORAPUCPRM)     /    CORMA)    670,675,675 
67C    IP     (EAcF-.2)     675,671,671 
«-7i    RUMP    ■    2.*BU¥P 
675    CALL    ABMAX(NFT£ ,CORAP , CORMA, ICPRM) 

IP  (KR(17))  680,680,685 
fPC  ic  (Koriq))  a?0,705,6^0 
6*5    CHNTTNJIF 

KR(  17)     =    K9(  17)     -    1 

««CC  comtinuf 

6l95    FORM/»T(2X23H0EBUG    CORRECTIONS    «=L,NL    ) 

WPITF(  KnijTT*,9c) 

WRTTF(KCL'Tf  i°0)  (FLE(  J)  ,  J=l  ,MAT1  I) 

WP!^P(KPUT,1?C)  (0VNL(  J I  ,J  =  l  ,MAT1J) 
705    CONTINUE 

CORRECT    PPI^APY    VARI^BL^S 

OUM=. 1 /RUM0 

EA«F=AvtNl(3.*F*SF.l.  ,OUW/ARS(CPOMA)  ) 
710    IF    (KR(13))     72C,720,715 
715    OHM    =    KP( 13) 

EA*C    =    fiviNKOUM    /    10..PASE) 
^20    IF    (EASE    -    1.0)     725,740,^40 
72*-    00    730    I  =  l,M*TH 
73C    FLp(  I  )     =    FLEU  )      *    EASE 

HP    7?5    I=1,MAT1J 
735    OVNL(I)     =    DVKL(I)     *    FACE 
^40    CONTINUE 
">50    00    765    J  =  2,NETA 

c(  1,J)     =    F(i  ,J)     +    PLF( J    -    1) 
7  ac    PPM  TlMtje 

M     =     jtlpT/ 

00    790    J=1,NETA 

M     =     M     +      1 

F( 2, J)    =    F(2,J)     ♦    FLF(M    -    l) 

78C    CONTINUE 

On    785    J=2,NFT£ 

M     =      M     +      1 

n*5    P(3,J    -    1)     =    F(3,J    -    1)     ♦    PLE(M    -    1) 
3CC    Pf 1,1 )=F( l,l)  +  PVNLU ) 

P(  3,NFTA  )=P(  3,NFTA)+0VNL(2) 

On    035    j=iTNETA 

F(4,  J)=F(4,J)+DVNLU+2) 

q-»5     rOMTTNUF 

IP     (  ITS    -    ^Q)     OCC,840,850 
9*0    Tc    (T7-77    _    777)     845,P50,R45 
045    T777    =    777 

ITc  =  3C 
«5C  CONTINUE 

RETURN 

FNO 

SUBROUTINE  I  NONE 

OIMPNSTPN  RRP(5) 

C0MM0N/ERPCn*/FLE<73)  ,FLA(^3) , F LEMf ELM C 3 ) , EL MM, IFLM,NE 
lLM,ILMMfOFL(73)  ,FNLE(28) ,ENL(2B  > ,CNLE¥, ENLM( 3  ) , FNL^M, I 
2FNLM,NFNLM,INLMM,peNL(28) , 0RNL(6) 

CPMM0N/FTAC0M/ETM25)  ,0ETA(24)  ,OSQ(2M  ,  OCU  (  2<*  ) ,  8 1  (  24  ) , 
18  2(24),LAP< 1C0)1BA1(73,2P) 

COMMCN/HTSC0M/TZ.T1 ,T2,TF(25) ,C 1 ,C2 ,C3 , C< ,BFT A, ZM( 5, 25 
l),XT( c0)  ,HP(25,5),HUP,HHUF,0LX2,C3M(50f ,*ETAM($0> 

COMMON/INTCPM/KR(2  0)  , KI N ,K0UT , MAT1 1 . MAT  I  J . NFT A, I ,  IS,NS 
l,IT,NTTME,NAMfMLE0,NNLP0,NONL1ITS,CASE(15),B(3i,MWFfN0 
2N,ITEMfMITFM,KR17  tNRx.N3T2,I0ENT,UMFS,FRF0tEPS 

COMMPN/MnNC0M/AM(2  3,2  8) ,0VNL(2« ) , C5 , C6 , C7 , C8 , CO , C 10, CI 
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ll,C12 

C0MM0N/PRPC0M/PR,SC,XM(5) 

COMMON/VARCOM/F(5,25) ,FTM1 (5  ,25 )  ,FTM2 ( 5 , 25 ) 
C  EVALUATF    GROUPINGS    WHICH    CONTRIBUTE    TO    (1-1)     PORTION 

C  OP    COEFFS    VARIABLES    WITH    DIMENSION     (NETA-1) 

T*=( I-2KC1, 4002,401 
4  0C2    XM(5)=C. 
C  EVALUATE    XM.    WHICH    CONTRIBUTES    TO    ERRORS    AND    TO    COEFFS 

C  AT    (  I  )    *ND    (  1-1) 

401    CALL    TAYLOR     (OF TA ( 1-1 )  , F (2  , 1-1 )  ,F( 2 , I ) , X M) 
C  FVAL    PORTION    OF     NLE    DEPENDENT    ON    XMt...    AND    GROUPINGS 

C  EVAL    AT    (1-1). 

40  3    C72=F( 2,1 )*X*(l)+F(3,T)*XM(2)*F(4v I )*XM(3 )+F(4, I-1)*XM 
1(4) 
XM( 5)=XM(5)+C72 

C73=F(  2,T)*7M(l,I-l)+F(3,I)*7.M(2,I-l)*FUfI)*ZM(3,  1-1) 
H-c(4,  I-1)**MU.I-1 ) 

CNLF( I  )=-(-C7-Cll+DETAU-l)/2.*(8ETA+C^*C'5)  +  .5*C2*C72- 
12.*C73) 
RP<M 1)=FNLE( I) 
C  FVAL    °ORTION    OF    ORIG    COEFFS    OF    AM    DEPENDENT    UPON    PAPAM 

C  FVAL    4T    U-l). 

*C5   AM(If I+1)-C2*XM(4)-2.*ZM(4,I-1I 
4C7    I»=(  1-2)    418,417,418 
^17    AM( 2,1)=-CC 
C  FVAL    1-1    PORTION    OF    COEFS,    AND    FORM    CONTRIBUTIONS    TO 

C  AM    AND    FRRORS. 

419    L°I=NETA-2+I 

DUM1=-C1C 
*21    CALL    AMSET    ( DUM1 , L PI  ,2  ,RPR ) 
411    LRI=LRI+NET^ 

DUM1=-1. 
A23    CALL    AMSFT    ( DUM1 ,L PI  ,3  ,RPR ) 
TP( I-2)*24,42*,424 

424  CALL    AVS^T    (  -CQ  ,  1-2  ,A  ,RRR) 

425  CONTINUE 
RRR(5)=FNLE(  I) 

oq9q    PnRMAT(2X,I5,2XlR5E12.4) 

IP(*R(17))     9901, 9902,  9901 
9901    CONTINUE 

WR!TF(KOUT,9G01) 
°001    F0RMAT(34H       DEBUG    FNLE    FRO*    IMCNE    PLUS    COEFS     ) 

WRITF(K0UT,Q999) I ,RRR 
9°02    CONTINUE 

pcTypNj 

END 

SUBROUTINE    A*SFT    ( CAM, LPI , NN, RR R ) 

DIMENSION    RRR(5) 

COMMON /EPPC0V/FLF(73)  ,ELA(73)  ,F LEM, ELM ( 3 ) , 5LMP,  IFLM,NE 
1LM,IL"M,PFL(73) ,FNLE(2  8) ,ENL(2*) ,PNLEM, ENLM( 3 ) , ENLMM, I 
2FNLM,NENLMt!NLMM    DFNL(28) ,DRNL(6) 

CnMM0N/FTACOM/ETA(25) ,DETA ( 24 > , D$Q( 2M , DCU ( 24 ) , Bl ( 24 ) , 
IB 2 (2*>  ,LAR(1C0) ,BA1( 73,28) 

C0MMnN/INTC0M/KR(20)  , KI N  ,KOUT , MAT1 I ,  MATU,NETAf I, IS,NS 
l,IT,NTIMr  ,NAV,NLEC,NNLEQ,NRNL,ITS,CA<>E(15),B(8),MWE,NO 
2N,!TEMfNITEM,KR17,NBT,NBT2  , I  DENT ,UMF^ , FRFQ, EPS 

C0MM0N/N0NCOV/AM( 2  8,28) ,0VNL (2P ) ,C5 , C6 ,C7, C8 , C9, C 10, 
1C11.C12 

FNLE( T  )=FNLE ( I ) -CAM*FLF ( LP  I ) 

IP(NN-1)2,1,1 

1  RRP(NN)=^NLE(I ) 

2  no   3   j=1,matij 

3  AM(  I, J)=AM(I  ,J)-CAM*8A1 (LP I, J) 

RcT!jpN 
END 

SUBROUTINE  ICNIY 
DIMENSION  RRRC4) 
C0MMQN/ERRC0M/FLE(73)  ,ELA(73)  ,FLEM, ELM ( 3 ) , ELMN, IFLM,NE 
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1LM,ILMM,PFL(73)  ,FNLE(29)  ,ENL(  28  )  ,FNLEM,  EMM(  3  ) ,  ENLMM,  I 
2FNLFM,NFNLM,INLMM,DFNL(2R) ,DRNL(6> 
COMMON /^T&CGV/FTA (25)  ,0ETA(24) ,DSQ(24) , DCU ( 24 ) , Bl ( 24 ) t 
1R2(  24)  ,1 AR(IGO)  ,BA1(73,28) 

COMMON/HISCOM/tj  n  ,T2  ,TF (2* ) ,C 1 ,C2 ,C3 ,C4, RET  A, ZM( 5, 25 
1),XH 50), HF( 25, 5) ,HUC ,HHUE ,0LX2 ,C3M ( 50 ) .PET AM (50) 

CPMMON/INTCnv/K«(2  0)  ,KINfKrtlT,MATlT.MATlJ.NFTAf  I,  IS,  NS 
l,r.MT]vr  ,N4M,NlFQ,NNLEQ,NRNL,ITS,CASE(  15  ) ,  B  (  8)  ,  MW  =  ,  NO 
2N,TTcwtNTTFMfKR17tNBT,NBT2  ,  I  PENT  ,  ,UMFf,  , FREQ,FPS 

COMMON /NPNCOM/AM( 2  8,2  8) ,DVNL(28) ,C5 , C6 , C7 , C8 , C^, C 10, 
ICU,C12 
C0MM0N/PPPC0M/PP,SC,XM(5) 

rnMM0N/vaRC0M/F(e,25) ,FTM1 (5  ,25 ) ,FTM2(5,25) 
C  ADn    rpNTOXRijTlONS    OF     I    T0    NONLINEAR    ERRORS 

C  EVALUATF    GROUPINGS    WHICH    ARE    USED    ONLY    AT     I     (NOT    AT 

C  1-1). 

FN1  F(  I  )=FNLE(I )-(C7  +  DETA(I-l) /2 .*< BFTa  +  c* *C5 ) +C 1 1 ) 
RRR(  i  )=FMLF( I) 
C  FVAL    °OPTTON    OF    PRIG    CPEFF^    OF     S«    Pfcf-njofNT    UPON    PARAM 

C  F  VA  L    A  t    t  # 

40  3    fiM(  I  T  I+2)=*M(I  ,I  +  2)+C2*XM(  3) -2.  *ZM< 3,1-1) 
C  AOP    CONTPIRUTICN^    OF     I     OEPcNOENT    CPEFFS    TO    AM    \ND 

C  NONLTNF^R    «=RRO°$. 

IF(  I-\'c^i)43<:,^38,430 

43  8    £M(  I,  2)=AMM  ,2)+l.+C2*XM(2)-2.*7M(2,I-l) 
43=    CALL    AMSET    (  C<*  ,1 -1  ,2  ,RPR) 

LP*=NFTA-H-I 

P'JM1=C  1C+C2*XM(  l)-2.  *ZM(  1,1-1) 
444    CALL    AMSFT    ( DUM1 , L PI  , 3  ,  RRR ) 
418    IF( I-NFTA )4KCf^Cc,459 
A«50    LP!  =  2*NFTA-l-ff 

DUM1=1.+C2*XM(2)-2.*?M(2,I-1) 
461    CUL    ^?H    (0UM1  ,LPT  ,4, RRR) 
405    CPNTINUF 
990c    F^RM^T( 2X,I5,2X1P4E12.4) 

I<=(KR(17»)     9  =  01  ,F902,<5^01 

Q«01    CONTINUE 

WRTTFftcnyT^cQ!) 

OQCl    FnpM*T(3^H       OEPUG    CNLE    FROM    IONLY    PLUS    CPEFS     ) 

V/DITp      fROU""",  Cqoq)      T,RRR 
9902    C0NTIM|)F 

RP^IJPN 

FNP 

SUBROUTINE    OL'TDIJT 

RC»1_  MU 

OIMENSION  CM(5)  ,U(25)  ,Y(25) 

COMMCN/crGCOM/rEM(^c) ,PTE(^C) ,OUES ,UE( 50 ,50 ) , °F(50,C0) 
l,RH0,MU,TtOSIP(r0)  ,IPSIP 

COMMON  /FTACOM/ETA( 25)  ,DETA ( 24) , OSQ( 24 > , DCU ( 24  ) , Bl ( 24  ) , 
18  2  <  24)  ,l?R(  100)  ,BA  1(73,2°) 

COMMON /HI SC0M/TZ,T1 ,T2,TF(25)  ,C1,C2,C3,C- .BETA, ?M(5,25 
1)  ,XT(  5C)  ,HF(25,5)  ,HUF,HHUE,PLX2  ,C3  M  (*Q  )  ,  RET  AM  (  50  ) 

COMMON /INTCOM/KR  (2  0)  ,  K  I  N  ,KOi  JT  ,  M  ATI  I  ,  MA^  l j , NET  A,  I,  IS,NS 
1,  IT,  NTT  Mr  ,NAM,NIfQ,NNLFQ,NRNL.TtS,CASE(  lc  )  ,  B  (  «  )  ,  MWF,  NO 
2V,ITFW,KIT'M,KPl7,NBT,NBT2,nENT,UMFS,FRco,FRS 

CnMM0N/N0NCnv/AM(2  8,2  8) ,DVNL (2° ) ,C5 , C6 , C7 ,CP , C9 , CIO, 
1C11,C12 

CPMMON/FRMCOV/TIME (CC) ,PRE ( 50 ) , PTET , S ( *0 ) , ROK A" ( 50 ) , RN 
10^F,VKACfN0TSC ,  I  CISC (50)  ,NSO(10) , MS  0(10) , I DPF,RA0NO, CO 
2N^,0ELT, CYCLE 

COMMON /PR  PC C*/PR,S0,XM(c) 

COMMON /TFMC0V/S°DUM2  (10)  ,DFR(5C)  ,0UMM1  (  2*  ) ,  S  LOPF  (  25  ) ,  R 
1F0HM(  2C)  ,S01JM1  (c0)  ,SPIIM2(50)  ,FWDUM(50)  ,  X  I  CON  (50  )  ,  FWCO 
2N(  c0)  ,PljnS(5C) 

COMMON/VSRCOM/F(5,25) ,FTM1(5,25) ,FTM2(5,25) 

COMMON /WALCOM/FW( 50) ,RHOVW(50) , TFW , I RHOVW , JRHGVW 
1    FOPMATC /°X65HFOGE  STATIC  XI 

1  FLUX    NORM  ROKAP    /6X57HVEL0C ITY  PRESSURE 

2(LB/  BFTA  PARAMETER    /6X*-7H(  FT/SCC  ) 

3(*TM)  SEC)**2  (FT) 
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WALL 


A    /3X6( 2X1PE1C.3) /) 

2  FORMAL /10X54HWALL    FR.  MOM    TRANS 

1  M*SS    FLUXES    /11X56HCOEFF,  COEFF,  SHEAR 

2  WFXH  TOTAL    /13X53HCF  RH0*UE*CF/2 

3  REM  r,AS    /3  =  X31H(LB/SQ    FT) 
A       (LP/SQ    FT    SEC)     /6X5 (3XIPE10 .3 ) / ) 

3  FORMA T( /10X58HREYMOLDS  STROUHAL  MOMENTUM  D 
IIS^LACF.  SHAPE  /UX58HNUMBER  NUMBER  THIC 
2KNFSS,          THICKNESS,            FACTOR,    /10X60H 

3  THETA  DELSTAP  DELSTAR/    /38X30HCF 

4T)  (FT)  THETA    /6X5C3X1»E10.3)/I 

A    FORMAT(/31X17HNPDAL    INFORMATION    // 3X50HDISTANCE 

1  STREAM  VELOCITY  SHEAR    /2X75HFROM   WALL 

2  FTA  ^UNCTION  RATIO  FUNCTION  FPPP 

3  VELOCITY    /5X72H(FT)  (F)  (FO=(j/UE 
A)             (FPO)                                             (FT/SEC)     /) 

5    FORMAT(1X7(1X1PE10.3) ) 
*    chdmat{  1PPF1C.3) 

DATA     IBUNK/2H       / 

IF    (KR(ll)-l)    300,300,301 
3C0    Ic    (KR(C))3Q2,3C2,301 
3C2    "H0VW1  IS)=C1*F(  1»1>+HF(1  ,5) 
3C1    CPP=-C3M(  IS)*MU 

0UM1=-1./C3M(IS) 

0UM2=RH0VW(I S) /C3M(I  S) 

IF    (UE( IS, ITEM) -l. ) 301 2, 301 1,30 12 

3011  UE( IS,ITFM)=0. 

3012  CONTINUE 
3051  CONTINUE 

DEP( 1)=C.C 

OEP( 2)=C.C 

VMECH=-rUM2 

DUM^=VMECH*1C0. 

IF    (0L'WA-0UM2)     1°01, 190,190 
19C1    VMPCH=0. 
190    IE(A8S(BETA)-.00C1)     303,30A,30< 
3C3    BETA=C. 
3CA    Y( 1)=0. 

DO    182    I=2,NETA 
182    Y( I )=Y( I-1)+C89*0ETA(I-1) 

SHEAR =MU*UE( IS, ITEM) /CP9*F(3,1) /32.1^A05 

WR!TE(KOUTfl)UE(IS,ITEM)  ,PE(IS,ITEM) ,X  I  dS  ) ,  BET  A  ,  DUM  1, 
lROKAOf  TS) 
212    ^UM1=PH0*UE(  IS,ITEM)  /*li 

C«r  =  MU/CP9*F(3,l) 

0UMA=1. 

21A   R»=S=nuvi*S(I  S) 

WCC=CF*2. /RHG/UECI S,!TEM) 

DUM3=CQ9*(FC l,MFTA)-F(l ,1) ) 

OELST=Y(NETA)-DUM3 

REDFLS=0UM1*DELST 

THMOM=OUv3-C89*XM(5) 

RCTHMn=n»JMl*THMOM 

SHAPF=OELST/THwry 

c;TR  =  PREO/UE(  IS, ITEM) 

OUMi=nijMl*s(TS) 

STR=*TR*S< I S) 

WR!TF(KCUT,2)WFC  ,CF  ,  SHEAR,  VMECH.DIJM2 

WPTTF(  KEUT,3)DUM1  ,  <:TP,THM0W,  DEL*  T,  SHAPE 

WP!TP(KOUT,A) 


00     183    !=1,NETA 
COMPUTE    TRUE    VALUES    OF 
REF1=E( 1,1 


F(I  ,J)     AND    ETA 


183 

1930 

1<34 


PEF2=F(2,I 
REF3=F(3,I 

RCF4=F(4,I 

PETA=ETAU 

U(  !)=PEF2*UE(IS,ITEM) 

=TAN=RETA /CFTA(NETA) ) 

WRITE (KCUT.5) Yd )  ,RETA,REF1,REF2,REF3,RFFA,U( I) 

IE    { IDENT-IRLANK)     19A,309,19A 

IE(M0N)    3C9,308,309 
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3CP 

3C9 

31C 
311 


r)IIM3  =  RH0*UE(  IS,  ITEM) 

DUW4=-1.  /C3 

CONTINUF 

IPf  KP(  1M  MIC, 311. 310 

WRTTF(Krt)T,^)(  (F(t  ,.J)  f  J  =  1,W=TA)  ,1  =  1,4) 

RCTIJPM 


95 


LIST  OF  REFERENCES 


1.  Schlichting,  H. ,  Boundary  Layer  Theory ,  4th  Ed., 

p.  107-265,  McGraw-Hill,  1960. 

2.  Lighthill,  M.  J.,  "Response  of  Laminar  Skin  Friction 

and  Heat  Transfer  to  Fluctuations  in  Stream  Velocity," 
Proceedings  of  the  Royal  Society  of  London,  A  224, 
June  1954. 

3.  Glauert,  M.  B. ,  "The  Laminar  Boundary  Layer  on  Oscil- 

lating Plates  and  Cylinders,  "Journal  of  Fluid 
Mechanics ,  v.  1,  Part  1,  May  1956. 

4.  Carrier,  G.  F.  and  DiPrima,  R.  C,  "On  the  Unsteady 

Motion  of  a  Viscous  Fluid  Past  a  Semi-Infinite  Flat 
Plate,"  Journal  of  Mathematics  and  Physics,  January 
1957. 

5.  Nickerson,  R.  J.,  The  Effect  of  Freestream  Oscillations 

on  the  Laminar  Boundary  Layer  on  a  Flat  Plate,  D.Sc. 
Thesis,  Massachusetts  Institute  of  Technology,  1957. 

6.  Hill,  P.  G.  and  Stenning,  A.  H.,  "Laminar  Boundary 

Layers  in  Oscillatory  Flow,"  Transactions  of  the  ASME: 
Journal  of  Basic  Engineering,  A  265,  September  1960. 

7.  Koob ,  S.  J.  and  Abbott,  D.  E.,  "Investigation  of  a 

Method  for  the  General  Analysis  of  Time  Dependent 
Two-Dimensional  Laminar  Boundary  Layers,"  Transactions 
of  the  ASME:  Journal  of  Basic  Engineering,  NO.  4,  1968. 

8.  Los  Alamos  Scientific  Laboratory  Report  LA-2910,  A 

Method  for  Computing  Nonsteady,  Incompressible, 
Viscous  Fluid  Flows,  by  J.  E.  Fromm,  May  196  3. 

9.  Chorin,  A.  J.,  "A  Numerical  Method  for  Solving  Incom- 

pressible Viscous  Flow  Problems,"  Journal  of  Compu- 
tational Physics,  v.  2,  1967. 

10.  Smith,  A.  M.  O.  and  Clutter,  D.  W. ,  "Solution  of  the 

Incompressible  Laminar  Boundary  Layer  Equations," 
AIAA  Journal,  v.  1,  No.  9,  September  1963. 

11.  Aerospace  Research  Laboratories  Report  ARL  67-0050, 

Finite  Difference  Solution  of  the  Laminar  Boundary 
Layer  Equations  in  the  Von  Mises  Variables  with 
Applications ,  by  G.  Kleinstein,  March  1967. 


96 


12.  Aerospace  Research  Laboratories  Report  ARL  67-009  8, 

Operation  Instructions  and  Description  of  Computer 
Program  for  Finite  Difference  Solution  of  the  Laminar 
Compressible  Boundary  Layer  Equations  in  the  Von 
Mises  Variables,  by  F.  Nabi ,  May  196  7. 

13.  Aerospace  Research  Laboratories  Report  ARL  67-0045, 

Solution  of  the  Laminar  Boundary  Layer  Equations 
with  Mass  Transfer  by  the  GKD  Multimoment  Method, 
by  H.  E.  Bethel,  March  1967. 

14.  Aerospace  Research  Laboratories  Report  ARL  67-0052, 

An  Approximate  Solution  of  Von  Mises '  Boundary  Layer 
Equation ,  by  A.  Reggiori ,  March  1967. 

15.  Aerospace  Research  Laboratories  Report  ARL  67-0124, 

On  an  Alternative  Method  of  Numerically  Evaluating 
the  Linearized  Problem  of  an  Impulsively  Moving  Flat 
Plate  in  a  Rarefied  Gas  Flow,  by  K.  S.  Nagaraja  and 
J.  P.  Hudson,  June  1967 

16.  Aerotherm  Corporation  Report  No.  66-7,  Part  III, 

Nonsimilar  Solution  of  the  Multicomponent  Laminar 
Boundary  Layer  by  an  Integral  Matrix  Method,  by 
E.  P.  Bartlett  and  R.  M.  Kendall,  NASA  CR-1062, 
June  196  8. 

17.  Howarth,  L.,  "On  the  Solution  of  the  Laminar  Boundary 

Layer  Equations,"  Proceedings  of  the  Royal  Society 
of  London,  A  164,  January-February  1938. 

18.  Aerotherm  Corporation  Report  No.  6  8-42,  User's  Manual 

Boundary  Layer  Integral  Matrix  Procedure  (BLIMP) ,  by 
E.  P.  Bartlett  and  C.  Deblaye,  October  1968. 


97 


INITIAL  DISTRIBUTION  LIST 


1.  Defense  Documentation  Center 
Cameron  Station 
Alexandria,  Virginia  22314 

2.  Library,  Code  0212 
Naval  Postgraduate  School 
Monterey,  California  93940 

3.  Commander,  Naval  Air  Systems  Command 
Department  of  the  Navy 
Washington,  D.C.  20360 

4.  Dean  of  Research  Administration 
Naval  Postgraduate  School 
Monterey,  California  93940 

5.  Chairman,  Department  of  Aeronautics 
Naval  Postgraduate  School 
Monterey,  California  93940 

6.  Professor  James  A.  Miller 
Department  of  Aeronautics 
Naval  Postgraduate  School 
Monterey,  California  93940 

7.  LT  Barry  A.  Gastrock,  USN 
34  Reeceville  Road 
Coatesville,  Pennsylvania  19320 


No.  Copies 
20 


98 


UHCL&SSIEIED 


Security  Classification 


DOCUMENT  CONTROL  DATA  R&D 


{Security  classification  ol  tltla,   body  of  abstract  and  indexing  annotation  must  be  entered  when  the  overall  report  is  classified) 


originating   activity  (Corporate  author) 

Naval  Postgraduate  School 
Monterey,  California  93940 


2».   REPORT    SECURITY    C  L  A  SSI  F  I  C  A  T  IOC 

Unclassified 


26.  GROUP 


REPORT  TITLE 


NONSIMILAR  SOLUTION  OF  THE  LAMINAR  BOUNDARY  LAYER 
IN  AN  OSCILLATORY  FLOW  BY  AN  INTEGRAL  MATRIX  METHOD 


DESCRIPTIVE  NOTES  (Type  ot  report  and,  inclusive  dates) 

Aeronautical  Engineer's  Thesis ,  June  1969 


5     authoriS)  (First  name,  middle  initial,  laal  name) 


Barry  Allen  Gastrock,  Lieutenant,  United  States  Navy 


IEPORT    DATE 

June    1969 


7a.     TOTAL    NO.    OF    PAGES 

98 


7b.    NO.    OF    REFS 


18 


6a.    CONTRACT    OR    GRANT   NO. 


b.    PROJECT    NO. 


9a.    ORIGINATOR'S   REPORT   NUMBER(S) 


numbera  that  may  ba  at  signed 


DISTRIBUTION    STATEMENT 


Distribution  of  this  document  is  unlimited, 


II.    SUPPLEMENTARY    NOTES 


12.    SPONSORING   MILI  TAR  Y    ACTIVITY 


Naval  Postgraduate  School 
Monterey,  California  93940 


13.  ABSTRAC  T 


The  development  of  a  numerical  procedure  for  the 
treatment  of  nonsimilar,  unsteady,  laminar  boundary  layers 
is  presented.   The  solution  is  obtained  from  the  laminar, 
isothermal,  incompressible  boundary  layer  equations  employ- 
ing a  modification  of  the  integral  matrix  procedure  of 
Bartlett  and  Kendall.   Solutions  of  example  problems  are 
presented  for  steady  Blasius  and  Howarth  flows,  and  for 
oscillating  Blasius  flow.   Agreement  with  the  known 
classical  results  is  satisfactory  and  establishes  the 
general  feasibility  of  the  method.   Core  storage  require- 
ments of  130,000  bytes  allow  consideration  of  as  many  as 
25  nodal  points  across  the  boundary  layer,  50  time  incre- 
ments per  oscillation  cycle  and  50  streamwise  stations. 
Solution  of  oscillating  Blasius  flow  considering  8  nodal 
points  and  16  time  increments  requires  13.49  seconds  for 
one  streamwise  station  utilizing  IBM  360/67  time  sharing 
capabilities . 


DD 


'•"•  1473 

i  nov  es  I  "T  /  *J 
S/N    0101 -807-681 1 


(PAGE     1) 


99 


UNCLASSIFIED 

Security  Classification 


UNCLASSIFIED 


Security  Classification 


KEY     WORDS 


OLE  WT 


LINK     B 


ROLE  W  T 


BOUNDARY  LAYER 
LAMINAR  BOUNDARY  LAYERS 
BOUNDARY  LAYER  SOLUTION 
BOUNDARY  LAYER  COMPUTER  PROGRAM 
UNSTEADY  BOUNDARY  LAYER 


DD  ,F°1"..1473   back, 


101  -  SO  7- 68  2 ' 


100 


UNCLASSIFIED 


Security  Classification 


Ctaulora  == 
SHELF   BINDER 

■  Syroeuse,  N.  Y. 

■  Stockton,  Colif. 


thesG2357 

'J2SS  SOlU,IOn  of  the  lamin^  bound 


